home *** CD-ROM | disk | FTP | other *** search
/ Power-Programmierung / CD1.mdf / forth / compiler / fpc / source / complex.seq < prev    next >
Encoding:
Text File  |  1988-11-08  |  3.1 KB  |  106 lines
  1. \ Complex Floating Point Routines   by Mark Smiley
  2. \ Corrections and additions by R. L. Smith
  3. \ Complex numbers are entered as ordered floating point pairs
  4. \ with the real part second on the f.p. stack, and the
  5. \ imaginary part at the top.  For magnitude and angle, the
  6. \ magnitude is second on the stack and the angle is at the
  7. \ top of the stack.
  8.  
  9. ANEW COMPLEX
  10. FLOATS
  11. HEX
  12.  
  13. : ZDROP   ( Z: z -- )
  14.         FDROP FDROP ;
  15.  
  16. : ZDUP   ( Z: z -- z z )
  17.         FOVER FOVER ;
  18.  
  19. : ZSWAP   ( Z: z1 z2 -- z2 z1 )
  20.         FSWAP 3 FNSWAP FSWAP 2 FNSWAP ;
  21.  
  22. : ZNIP   ( Z: z1 z2 -- z1 )
  23.         FROT FNIP FROT FNIP ;
  24.  
  25. : ZNEGATE ( F: x y --  -x -y )  ( Z: z -- -z )
  26.         FNEGATE FSWAP  FNEGATE FSWAP  ;
  27.  
  28. : ZCONJUGATE   ( Z: z1 -- z2 )
  29.         FNEGATE ;
  30.  
  31. : ZJ*   ( Z: z1 -- z2 )  ( Multiply by j , the unit imaginary )
  32.         FNEGATE FSWAP ;
  33.  
  34. : Z+  ( Z: z1 z2 -- z3 )  ( adds 2 complex numbers )
  35.     FROT F+  F-ROT  F+  FSWAP  ;
  36.  
  37. : Z-  ( Z: z1 z2 -- z3 )  FROT FSWAP F- F-ROT F- FSWAP  ;
  38.  
  39. : Z*  ( Z: z1 z2 -- z3 )  ( Multiply 2 complex numbers )
  40.         FOVER 3 FPICK F* FOVER 5 FPICK F* F+ FPOP
  41.         FROT F* F-ROT F* F\- FPUSH ;
  42.  
  43. : Z/   ( Z: z1 z2 -- z3 )  ( Divide z1 by z2 )
  44. \ Based on my (R. L. Smith) algorithm in ACM
  45. ( Vol.6, No. 8 , 1962 )
  46.         FDUP0=
  47.         IF      FDROP FDUP0=
  48.                 IF      [ LAST @ NAME> ] LITERAL 1 FPERR
  49.                 ELSE    FPCOPY F/ FSWAP FPUSH F/ FSWAP
  50.                 THEN
  51.                 EXIT
  52.         THEN
  53.         FPSEXP 7FFF AND FOVER FDUP0=
  54.         IF      FDROP FNIP FPCOPY F/
  55.                 FSWAP FPUSH F/ FNEGATE DROP EXIT
  56.         THEN
  57.         FPSEXP 7FFF AND < 0=  ( Comp. magnitudes of exponents )
  58.         IF       ( Real part of denom. < imag. part of denom. )
  59.                 FDROP ZJ* ZSWAP ZJ* ZSWAP FOVER
  60.         THEN    ( Now imaginary part <= real part )
  61.         FOVER FSWAP F/                         ( F: A B C D R )
  62.         FPCOPY F* F+                           ( F: A B C+D*R )
  63.         F-ROT FDUP 2 FPICK FPUSH FPCOPY F* F-
  64.         F-ROT FPUSH F* F+ FROT FPCOPY F/    ( F: den qr B-A*R )
  65.         FSWAP FPUSH F/ ;
  66.  
  67. : ZVARIABLE   CREATE  16 ALLOT  DOES>  ;
  68.  
  69. : Z!  ( F: r1 r2 -- ; addr -- )  FSWAP DUP F!  8 +  F!  ;
  70.  
  71. : Z@  ( F: -- r1 r2 ; addr -- )  DUP  F@  8 +  F@   ;
  72.  
  73. : ZMAG   ( F: r1 r2 -- r3 )
  74.         FABS FDUP0=
  75.         IF      FDROP FABS EXIT  THEN
  76.         FPSEXP FSWAP FABS FDUP0=
  77.         IF      DROP FDROP EXIT  THEN
  78.         FPSEXP <
  79.         IF      FSWAP  THEN
  80.         FOVER F/ FDUP F* F1.0 F+ FSQRT F* ;
  81.  
  82. : TOPOLAR   ( F: x y -- r theta )
  83.         F2DUP FATAN FPOP ZMAG FPUSH ;
  84.  
  85. : TOCART   ( F: r theta -- x y )
  86.         F2DUP FSIN F* FPOP FCOS F* FPUSH ;
  87.  
  88. \ Complex Square Root Algorithm
  89.  
  90. : F^2  ( F: x -- x^2 )   FDUP  F*  ;
  91. : F2/  ( F: x -- x/2 )   FPOP 1- FPUSH  ;
  92. : F2*  ( F: r1 -- r2 )   FPOP 1+ FPUSH  ;
  93.  
  94. : ZSQRT  ( F: x y -- x1 y1 )
  95.         FDUP0< FOVER FDUP0< FOVER ZMAG
  96.         FROT FABS F+ F2/ FSQRT FSWAP FABS F2/
  97.         FOVER F/
  98.         IF      FSWAP   THEN
  99.         IF      FNEGATE THEN ;
  100.  
  101. : -ZSQRT ( F: x y -- x1 y1 )   ZSQRT  ZNEGATE  ;
  102.  
  103. DECIMAL
  104.  
  105.  
  106.