home *** CD-ROM | disk | FTP | other *** search
/ Frozen Fish 1: Amiga / FrozenFish-Apr94.iso / bbs / oct93 / develop / umbscheme.lha / UMBScheme / src / complex.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-08-04  |  13.8 KB  |  629 lines
  1. /* complex.c -- UMB Scheme, complex number package.
  2.  
  3. UMB Scheme Interpreter                  $Revision: 2.5 $
  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**2 onto stack */
  411.         Complex_Sqrt();             /* compute sqrt (1-z**2) */
  412.         Pop(1);                /* remove 1-z**2 */
  413.         Push( Value_Register );     /* push sqrt(1-z**2) */
  414.         Make_Complex_Number( -B, A);    /* make iz */
  415.         Push( Value_Register );     /* push iz */
  416.         Complex_Add();            /* compute iz + sqrt(1-z**2) */
  417.         Pop(2);                /* remove iz and sqrt(1-z**2) */
  418.         Push( Value_Register );     /* push iz + sqrt(1-z**2) */
  419.         Complex_Log();            /* compute log(iz + sqrt(1-z**2)) */
  420.         Pop(1);                /* remove iz + sqrt(1-z**2) */
  421.         Push( Value_Register );        /* push log(iz + sqrt(1-z**2)) */
  422.         Make_Complex_Number( D, -C);    /* compute -ilog(iz + sqrt(1-z**2)) */
  423.         Pop(1);                /* remove log(iz + sqrt(1-z**2)) */
  424.  
  425. }
  426.  
  427.  
  428.  
  429.  
  430. Public void Complex_Acos()
  431.     /* acosz = -ilog(z + i*sqrt(1-z**2)) */
  432. {
  433.         Make_Complex_Number(1.0-(C*C)+(D*D), -(2.0*(C*D)) ); /* make 1-z**2 */
  434.         Push( Value_Register );     /* push 1-z**2 */
  435.         Complex_Sqrt();             /* compute sqrt (1-z**2) */
  436.         Pop(1);                /* remove 1-z*2 */
  437.         Push( Value_Register );        /* push sqrt(1-z**2) */
  438.         Make_Complex_Number( 0.0,1.0);  /* make i */
  439.         Push( Value_Register );     /* push i */
  440.         Complex_Multiply();        /* compute i*sqrt(1-z**2) */
  441.         Pop(2);                /* remove i and sqrt(1-z**2) */
  442.         Push( Value_Register );     /* push i*sqrt(1-z**2) */
  443.         Complex_Add();            /* compute z + i*sqrt(1-z**2) */
  444.         Pop(1);                /* remove i*sqrt(1-z**2) */
  445.         Push( Value_Register );        /* push z + i*sqrt(1-z**2) */
  446.         Complex_Log();            /* compute log(z + i*sqrt(1-z**2)) */
  447.         Pop(1);                /* remove z + i*sqrt(1-z**2) */
  448.         Push( Value_Register );        /* push log(z + i*sqrt(1-z**2)) */
  449.         Make_Complex_Number( D, -C);    /* compute -ilog(z + i*sqrt(1-z**2)) */
  450.         Pop(1);                /* remove log(z + i*sqrt(1-z**2)) */
  451.  
  452. }
  453.  
  454.  
  455.  
  456. Public void Complex_Atan()
  457.     /* atanz = -ilog((1+iz)*sqrt(1/(1+z**2))) */
  458. {
  459.         Make_Complex_Number(1.0,0.0);    /* make 1 */
  460.         Push(Value_Register);        /* push 1 */
  461.         Make_Complex_Number(1.0+(A*A)-(B*B), 2.0*(A*B)); /* make 1+z**2 */
  462.         Push( Value_Register );     /* push 1+z**2 */
  463.         Complex_Divide();        /* compute 1 / (1+z**2) */
  464.         Pop(2);                /* remove 1 and 1+z**2 */
  465.         Push(Value_Register);        /* push 1/(1+z**2) */
  466.         Complex_Sqrt();             /* compute sqrt(1/(1+z**2)) */
  467.         Pop(1);                /* remove 1/(1+z**2) */
  468.         Push( Value_Register );        /* push sqrt(1/(1+z**2)) */
  469.         Make_Complex_Number( 1.0 - B, A); /* make 1+iz */
  470.         Push( Value_Register );        /* push 1+iz */
  471.     Complex_Multiply();         /* compute (1+iz)(sqrt(1/(1+z**2))) */
  472.         Pop(2);             /* remove 1+iz and sqrt(1/(1+z**2)) */
  473.         Push( Value_Register );     /* push (1+iz)(sqrt(1/(1+z**2))) */
  474.         Complex_Log();         /* compute log((1+iz)(sqrt(1/(1+z**2))) */
  475.         Pop(1);                /* remove (1+iz)(sqrt(1/(1+z**2))) */
  476.         Push( Value_Register );     /* push log((1+iz)(sqrt(1/(1+z**2))) */
  477.         Make_Complex_Number( D, -C);/* compute -ilog((1+iz)(sqrt(1/(1+z**2))) */
  478.         Pop(1);            /* remove log((1+iz)(sqrt(1/(1+z**2))) */
  479.  
  480. }
  481.  
  482.  
  483.  
  484. Public void Complex_Atan2()
  485. {
  486.     Error("No Atan2 Operation on Complex!");
  487. }
  488.  
  489.  
  490.  
  491.  
  492. Public void Complex_Exact_To_Inexact()
  493. {
  494.     Value_Register = Top( 1 );
  495. }
  496.  
  497.  
  498.  
  499. Public void Complex_Inexact_To_Exact()
  500. {
  501.         Error("inexact->to->exact is not applicable to complex numbers");
  502. }
  503.  
  504.  
  505.  
  506.  
  507. Public void Complex_To_String()
  508. {
  509.     Import    void    String_Append();
  510.     Integer    radix = Number_To_Integer( Top(1) );
  511.  
  512.     if ( A == 0.0 )
  513.     {
  514.         Make_Constant_String(     radix == 2 ? "#b" :
  515.                     radix == 8 ? "#o" :
  516.                     radix == 10 ? "" : "#x" );
  517.         Push( Value_Register );
  518.     }
  519.     else
  520.     {
  521.         Make_Real_Number( A );
  522.         Push( Value_Register );
  523.         Push( Top(2) );
  524.         Number_To_String(); Pop(2);
  525.         Push( Value_Register );
  526.     }
  527.  
  528.     if ( Get_Number_Complex_Imaginary_Part( Top(3) ) >= 0 )
  529.     {
  530.         Make_Constant_String( "+" );
  531.         Push( Value_Register );
  532.         String_Append(); Pop(2);
  533.         Push( Value_Register );
  534.     }
  535.  
  536.     Make_Real_Number( Get_Number_Complex_Imaginary_Part( Top(3) ) );
  537.     Push( Value_Register );
  538.     Push( Top(3) );
  539.     Number_To_String(); Pop(2);
  540.     Push( Value_Register );
  541.     if ( radix != 10 )
  542.     {
  543.         /* remove radix prefix from imaginary part */
  544.  
  545.         Make_Constant_String( Get_String_Value( Top(1) ) + 2 ); Pop(1);
  546.         Push( Value_Register );
  547.     }
  548.     String_Append(); Pop(2);
  549.     Push( Value_Register );
  550.  
  551.     Make_Constant_String( "i" );
  552.     Push( Value_Register );
  553.     String_Append(); Pop(2);
  554. }
  555.  
  556. Public void Complex_Make_Rectangular()
  557.         /* Make a Complex in rectangular form; z = x + yi */
  558. {
  559.         Tower_Position argtypeX = Get_Number_Tower_Position(Top(2));
  560.         Tower_Position argtypeY = Get_Number_Tower_Position(Top(1));
  561.  
  562.         if (argtypeX == COMPLEX_LEVEL || argtypeY == COMPLEX_LEVEL)
  563.         {
  564.                 Error("Bad arguments to Make-Rectangular.");
  565.         }
  566.  
  567. }
  568.  
  569.  
  570.  
  571. Public void Complex_Make_Polar()
  572. {
  573.         Tower_Position argtypeX = Get_Number_Tower_Position(Top(2));
  574.         Tower_Position argtypeY = Get_Number_Tower_Position(Top(1));
  575.  
  576.         if (argtypeX == COMPLEX_LEVEL || argtypeY == COMPLEX_LEVEL)
  577.         {
  578.                 Error("Bad arguments to Make-Rectangular.");
  579.         }
  580. }
  581.  
  582.  
  583.  
  584. Public void Complex_Real_Part()
  585.         /* Return the real part of a complex in rectangular form */
  586. {
  587.         Make_Real_Number( C );
  588. }
  589.  
  590.  
  591.  
  592.  
  593. Public void Complex_Imaginary_Part()
  594.         /* Return the imaginary part of a complex in rectangular form */
  595. {
  596.         Make_Real_Number( D );
  597. }
  598.  
  599.  
  600.  
  601.  
  602. Public void Complex_Magnitude()
  603.         /* Return the Magnitude of a complex, */
  604. {
  605.         Make_Real_Number( sqrt( (C*C) + (D*D) ) );
  606. }
  607.  
  608.  
  609.  
  610.  
  611. Public void Complex_Angle()
  612.         /* Return the angle of a complex */
  613. {
  614.     Double ang;
  615.     
  616.     ang = ( C == 0.0 && D > 0.0 ) ? HALF_PI :
  617.              ( ( C == 0.0 && D < 0.0 ) ? -HALF_PI : atan2(D,C) );
  618.  
  619.     if (Is_Complex_Zero()) Make_Real_Number( 0.0 );
  620.         else Make_Real_Number( ang );
  621.  
  622. }
  623.  
  624.  
  625. #undef A
  626. #undef B
  627. #undef C
  628. #undef D
  629.