home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 8 / FreshFishVol8-CD2.bin / bbs / dev / umbscheme-2.12.lha / UMBScheme / src / complex.c < prev    next >
C/C++ Source or Header  |  1993-11-29  |  14KB  |  629 lines
  1. /* complex.c -- UMB Scheme, complex number package.
  2.  
  3. UMB Scheme Interpreter                  $Revision: 2.12 $
  4. Copyright (C) 1988, 1991 William R Campbell
  5.  
  6. This program is free software; you can redistribute it and/or modify
  7. it under the terms of the GNU General Public License as published by
  8. the Free Software Foundation; either version 1, or (at your option)
  9. any later version.
  10.  
  11. This program is distributed in the hope that it will be useful,
  12. but WITHOUT ANY WARRANTY; without even the implied warranty of
  13. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  14. GNU General Public License for more details.
  15.  
  16. You should have received a copy of the GNU General Public License
  17. along with this program; if not, write to the Free Software
  18. Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  19.  
  20. UMB Scheme was written by Bill Campbell with help from Karl Berry,
  21. Barbara Dixey, Ira Gerstein, Mary Glaser, Kathy Hargreaves, Bill McCabe,
  22. Long Nguyen, Susan Quina, Jeyashree Sivasubram, Bela Sohoni and Thang Quoc Tran.
  23.  
  24. For additional information about UMB Scheme, contact the author:
  25.  
  26.     Bill Campbell
  27.     Department of Mathematics and Computer Science
  28.     University of Massachusetts at Boston
  29.     Harbor Campus
  30.     Boston, MA 02125
  31.  
  32.     Telephone: 617-287-6449        Internet: bill@cs.umb.edu
  33.  
  34. */
  35.  
  36.  
  37. #include "portable.h"
  38. #include "eval.h"
  39. #include "object.h"
  40. #include "architecture.h"
  41. #include "number.h"
  42. #include "fixnum.h"
  43. #include "bignum.h"
  44. #include "rational.h"
  45. #include "real.h"
  46. #include "complex.h"
  47. #include "steering.h"
  48. #include "io.h"
  49. #include <math.h>
  50. #include <errno.h>
  51.  
  52. #define A (Get_Number_Complex_Real_Part(Top(2)))
  53. #define B (Get_Number_Complex_Imaginary_Part(Top(2)))
  54. #define C (Get_Number_Complex_Real_Part(Top(1)))
  55. #define D (Get_Number_Complex_Imaginary_Part(Top(1)))
  56.  
  57. #define PI    3.14159265358979323846
  58. #define HALF_PI    1.57079632679489661923
  59.  
  60.  
  61. /* Predicates. */
  62.  
  63. Public Boolean Is_Complex_Zero()
  64. {
  65.     return ( (C == 0.0) && (D == 0.0) );
  66. }
  67.  
  68.  
  69. Public Boolean Is_Complex_Positive()
  70. {
  71.         Error("No Positive/Negative Concept about Complex!");
  72.         return FALSE;
  73. }
  74.  
  75.  
  76. Public Boolean Is_Complex_Negative()
  77. {
  78.         Error("No Positive/Negative Concept about Complex!");
  79.         return FALSE;
  80. }
  81.  
  82.  
  83.  
  84. Public Boolean Is_Complex_Odd()
  85. {
  86.         Error("No Odd/Even Concept about Complex!");
  87.         return FALSE;
  88. }
  89.  
  90.  
  91.  
  92. Public Boolean Is_Complex_Even()
  93. {
  94.         Error("No Odd/Even Concept about Complex!");
  95.         return FALSE;
  96. }
  97.  
  98.  
  99.  
  100. Public Boolean Is_Complex_Exact()
  101. {
  102.     return FALSE;    /* Complex numbers are always inexact */
  103. }
  104.  
  105.  
  106. Public Boolean Is_Complex_Inexact()
  107. {
  108.         return TRUE;    /* Complex numbers are always inexact */
  109. }
  110.  
  111. /* Comparisons. */
  112.  
  113.  
  114. Public Boolean Complex_Less_Than()
  115. {
  116.         Error("No Inequality about Complex!");
  117.     return FALSE;
  118. }
  119.  
  120.  
  121.  
  122. Public Boolean Complex_Greater_Than()
  123. {
  124.         Error("No Inequality about Complex!");
  125.     return FALSE;
  126. }
  127.  
  128.  
  129.  
  130. Public Boolean Complex_Equal()
  131. {
  132.     return ( (A == C) && (B == D) );
  133. }
  134.  
  135.  
  136. Public Boolean Complex_Less_Than_Or_Equal()
  137. {
  138.         Error("No Inequality about Complex!");
  139.     return FALSE;
  140. }
  141.  
  142.  
  143. Public Boolean Complex_Greater_Than_Or_Equal()
  144. {
  145.         Error("No Inequality about Complex!");
  146.         return FALSE;
  147. }
  148.  
  149.  
  150. /* Arithmetic. */
  151.  
  152.  
  153. Public void Complex_Add()
  154. {
  155.     Make_Complex_Number( A + C, B + D );
  156. }
  157.  
  158.  
  159.  
  160. Public void Complex_Subtract()
  161. {
  162.     Make_Complex_Number( A - C, B - D );
  163. }
  164.  
  165.  
  166.  
  167. Public void Complex_Multiply()
  168. {
  169.     Make_Complex_Number( ( (A*C) - (B*D) ), ( (A*D) + (B*C) ) );
  170. }
  171.  
  172.  
  173.  
  174. Public void Complex_Divide()
  175. {
  176.     if (Is_Complex_Zero())           /* C+Di = 0+0i, it's on top(1) */
  177.         Error("Complex division by ZERO !!!!");
  178.     else
  179.         Make_Complex_Number(     ( (A*C)+(B*D) ) / ( (C*C)+(D*D) ),
  180.                          ( (B*C)-(A*D) ) / ( (C*C)+(D*D) ) );
  181. }
  182.  
  183.  
  184.  
  185. Public void Complex_Quotient()
  186. {
  187.         Error("No Quotient Operation on Complex!");
  188. }
  189.  
  190.  
  191.  
  192. Public void Complex_Remainder()
  193. {
  194.         Error("No Remainder Operation on Complex!");
  195. }
  196.  
  197.  
  198.  
  199. Public void Complex_Modulo()
  200. {
  201.     Error("No Modulo Operation on Complex!");
  202. } 
  203.  
  204.  
  205.  
  206.  
  207. Public void Complex_Negate()
  208. {
  209.     Value_Register = Copy_Object(Top(1), Complex_Size);
  210.  
  211.     Get_Number_Complex_Real_Part(Value_Register) =
  212.         - Get_Number_Complex_Real_Part(Value_Register);
  213.  
  214.  
  215.     Get_Number_Complex_Imaginary_Part(Value_Register) =
  216.         - Get_Number_Complex_Imaginary_Part(Value_Register); 
  217.  
  218. }
  219.  
  220.  
  221.  
  222. Public void Complex_Abs()
  223. {
  224.         Make_Real_Number(sqrt( (C*C) + (D*D) ));
  225. }
  226.  
  227.  
  228.  
  229. Public void Complex_Numerator()
  230. {
  231.     Error("No Numerator Operation on Complex!");
  232. }
  233.  
  234.  
  235.  
  236. Public void Complex_Denominator()
  237. {
  238.     Error("No Denominator Operation on Complex!");
  239. }
  240.  
  241.  
  242.  
  243. Public void Complex_Rationalize()
  244. {
  245.     Error("No Rationalize Operation on Complex!");
  246. }
  247.  
  248.  
  249.  
  250. Public void Complex_Max()
  251. {
  252.         Error("No Max Operation on Complex!");
  253. }
  254.  
  255.  
  256.  
  257. Public void Complex_Min()
  258. {
  259.         Error("No Min Operation on Complex!");
  260. }
  261.  
  262.  
  263.  
  264. Public void Complex_GCD()
  265. {
  266.     Error("No GCD Operation on Complex!");
  267. }
  268.  
  269.  
  270.  
  271. Public void Complex_LCM()
  272. {
  273.     Error("No LCM Operation on Complex!");
  274. }
  275.  
  276.  
  277.  
  278. Public void Complex_Floor()
  279. {
  280.     Error("No Floor Operation on Complex!");
  281. }
  282.  
  283.  
  284.  
  285. Public void Complex_Ceiling()
  286. {
  287.     Error("No Ceiling Operation on Complex!");
  288. }
  289.  
  290.  
  291.  
  292. Public void Complex_Truncate()
  293. {
  294.     Error("No Truncate Operation on Complex!");
  295. }
  296.  
  297.  
  298.  
  299. Public void Complex_Round()
  300. {
  301.     Error("No Round Operation on Complex!");
  302. }
  303.  
  304.  
  305.  
  306.  
  307. Public void Complex_Sqrt()
  308.     /* sqrtz = sqrt(|z|)*( cos((angle+2kPi)/2) + isin((angle+2kPi)/2) ) */
  309. {
  310.         Double mag_sqrt = sqrt(sqrt( (C*C) + (D*D) ));
  311.     Double ang = ( C == 0.0 && D > 0.0 ) ? HALF_PI :
  312.              ( ( C == 0.0 && D < 0.0 ) ? -HALF_PI : atan2(D,C) );
  313.  
  314.     /* select a proper angle so that 
  315.                 the real part of the result is non-negative */
  316.  
  317.         if (cos(ang/2.0) >= 0.0) ang /= 2.0;
  318.         else ang = (ang/2.0) + PI;
  319.  
  320.         Make_Complex_Number(mag_sqrt * cos(ang), mag_sqrt * sin(ang));
  321. }
  322.  
  323.  
  324.  
  325. Public void Complex_Exp()
  326. {
  327.         Make_Complex_Number( exp(C) * cos(D) , exp(C) * sin(D) );
  328. }
  329.  
  330.  
  331.  
  332. Public void Complex_Log()
  333.     /* logz = log(magnitude(z)) + angle(z)i */
  334. {
  335.     Double ang = ( C == 0.0 && D > 0.0 ) ? HALF_PI :
  336.              ( ( C == 0.0 && D < 0.0 ) ? -HALF_PI : atan2(D,C) );
  337.  
  338.     if (Is_Complex_Zero()) Display_Error("Bad arg to Log : ", Top(1) );
  339.     else  Make_Complex_Number( log(sqrt(C*C + D*D)), ang ); 
  340. }
  341.  
  342.  
  343.  
  344. Public void Complex_Expt()
  345.         /* z1**z2 = e**z2log(z1)   ;  0**0 is defined to be 1 */
  346. {
  347.  
  348. #define Z1_and_Z2_are_Zeros (A == 0.0 && B == 0.0) && (C == 0.0 && D == 0.0)
  349. #define Z1_is_Zero  A == 0.0 && B == 0.0
  350.  
  351.     Promote( 2 , COMPLEX_LEVEL );
  352.     if ( Z1_and_Z2_are_Zeros )  Make_Complex_Number(1.0, 0.0);
  353.     else if ( Z1_is_Zero )    Display_Error("Bad arg to Expt:", Top(2));
  354.     else 
  355.     {
  356.             Push(Top(2));         /* push z1 */
  357.             Complex_Log();      /* compute logz1 */
  358.             Pop(1);            /* remove z1 */
  359.             Push(Value_Register);    /* push logz1 */
  360.             Complex_Multiply();     /* compute z2logz1 */
  361.             Pop(1);            /* remove logz1 */
  362.             Push(Value_Register);    /* push z2logz1 */
  363.             Complex_Exp();        /* compute e**z2logz1 */
  364.             Pop(1);            /* remove z2logz1 */
  365.          }
  366. #undef Z1_and_Z2_are_Zeros
  367. #undef Z1_is_Zero 
  368. }
  369.  
  370.  
  371.  
  372. Public void Complex_Sin()
  373.     /* sinz = (e**iz - e**-iz) / 2i */
  374. {
  375.         Make_Complex_Number( ((exp(-D) * sin(C)) - (exp(D) * sin(-C))) / 2.0 , 
  376.                  ((exp(D) * cos(-C)) - (exp(-D) * cos(C))) / 2.0 );
  377. }
  378.  
  379.  
  380.  
  381. Public void Complex_Cos()
  382.     /* cosz = (e**iz + e**-iz) / 2 */
  383. {
  384.         Make_Complex_Number(((exp(-D) * cos(C)) + (exp(D) * cos(-C))) / 2.0 ,
  385.                             ((exp(-D) * sin(C)) + (exp(D) * sin(-C))) / 2.0 );
  386. }
  387.  
  388.  
  389.  
  390. Public void Complex_Tan()
  391.     /* tanz = sinz / cosz */
  392. {
  393.         Complex_Sin();          /* compute sin(z) */
  394.         Push( Value_Register ); /* push sin(z) onto stack */
  395.         Push( Top(2) );         /* swap z upto top to compute cos(z) */
  396.         Complex_Cos();          /* compute cos(z) */
  397.         Pop(1);                 /* remove the temp z */
  398.         Push( Value_Register ); /* push cos(z) onto stack */
  399.         Complex_Divide();       /* compute tan(z) = sin(z) / cos(z) */
  400.         Pop(2);                 /* remove sin(z) and cos(z) off the stack */
  401. }
  402.  
  403.  
  404.  
  405. Public void Complex_Asin()
  406.     /* asinz = -ilog(iz + sqrt(1-z**2)) */
  407. {
  408.         Make_Complex_Number( 1.0 - (C*C) + (D*D), -(2.0*(C*D)) ); 
  409.                         /* Make 1-z**2 */
  410.     Push( Value_Register );     /* push 1-z**