home *** CD-ROM | disk | FTP | other *** search
/ OS/2 Shareware BBS: 10 Tools / 10-Tools.zip / mod201j.zip / modula2.exe / os2src / complexl.mod < prev    next >
Text File  |  1995-12-26  |  16KB  |  465 lines
  1. IMPLEMENTATION MODULE ComplexLib;
  2.  
  3. (**************************************************************************
  4.   OS/2 2.x  Modula-2 complex number library
  5.  
  6.             Aug 16, 1995 -- Created.
  7.             Dec 25, 1995 -- Updated: No more FSAVE/FRSTOR needed; Neuhoff
  8.  
  9.   Copyright (c) 1995 by Anthony Busigin. Permission is granted for
  10.   use of this software with the Modula-2 OS/2-compiler by Juergen Neuhoff.
  11.  
  12.   Note: This library is a work in progress. SHORTCOMPLEX functions will
  13.         be added in a future version of the library.
  14. **************************************************************************)
  15.  
  16. (*$XL+*)
  17. (*$F+  allow coprocessor instructions for INLINE assembler *)
  18.  
  19. IMPORT
  20.   InOut,
  21.   MathLib0,
  22.   RealInOut,
  23.   SYSTEM;
  24.  
  25. FROM SYSTEM IMPORT INLINE, ADR;
  26.  
  27.  
  28. PROCEDURE cadd( x, y : LONGCOMPLEX ) : LONGCOMPLEX;
  29. (* returns x+y *)
  30. BEGIN
  31.   INLINE
  32.   (
  33.   PUSH    EDI
  34.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  35.   FLD     x.r[ EBP ]                    (* x.r *)
  36.   FADD    y.r[ EBP ]                    (* y.r+x.r *)
  37.   FSTP    QWORD PTR [ EDI ]             (* cadd.r := ST *)
  38.   FLD     x.i[ EBP ]                    (* x.i *)
  39.   FADD    y.i[ EBP ]                    (* y.i+xi *)
  40.   FSTP    QWORD PTR [ EDI + 8 ]         (* cadd.i := ST *)
  41.   FWAIT
  42.   POP     EDI
  43.   );
  44. END cadd;
  45.  
  46.  
  47. PROCEDURE csub( x, y : LONGCOMPLEX ) : LONGCOMPLEX;
  48. (* returns x-y *)
  49. BEGIN
  50.   INLINE
  51.   (
  52.   PUSH    EDI
  53.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  54.   FLD     x.r[ EBP ]                    (* x.r *)
  55.   FSUB    y.r[ EBP ]                    (* y.r+x.r *)
  56.   FSTP    QWORD PTR [ EDI ]             (* cadd.r := ST *)
  57.   FLD     x.i[ EBP ]                    (* x.i *)
  58.   FSUB    y.i[ EBP ]                    (* y.i+xi *)
  59.   FSTP    QWORD PTR [ EDI + 8 ]         (* cadd.i := ST *)
  60.   FWAIT
  61.   POP     EDI
  62.   );
  63. END csub;
  64.  
  65.  
  66. PROCEDURE cmul( x, y : LONGCOMPLEX ) : LONGCOMPLEX;
  67. (* returns x*y *)
  68. BEGIN
  69.   INLINE
  70.   (
  71.   PUSH    EDI
  72.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  73.   FLD     x.r[ EBP ]                    (* r *)
  74.   FLD     x.i[ EBP ]                    (* i r *)
  75.   FLD     y.r[ EBP ]                    (* rr i r *)
  76.   FLD     y.i[ EBP ]                    (* ii rr i r *)
  77.   FLD     ST                            (* ii ii rr i r *)
  78.   FMUL    ST, ST(3)                     (* i*ii ii rr i r *)
  79.   FXCH    ST(1)                         (* ii i*ii rr i r *)
  80.   FMUL    ST, ST(4)                     (* r*ii i*ii rr i r *)
  81.   FXCH    ST(2)                         (* rr i*ii r*ii i r *)
  82.   FMUL    ST(3), ST                     (* rr i*ii r*ii i*rr r *)
  83.   FMULP   ST(4), ST                     (* i*ii r*ii i*rr r*rr *)
  84.   FSUBP   ST(3), ST                     (* r*ii i*rr (r*rr-i*ii) *)
  85.   FADDP   ST(1), ST                     (* (r*ii+i*rr) (r*rr-i*ii) *)
  86.   FSTP    QWORD PTR [ EDI + 8 ]         (* cmul.i := ST *)
  87.   FSTP    QWORD PTR [ EDI ]             (* cmul.r := ST *)
  88.   FWAIT
  89.   POP     EDI
  90.   );
  91. END cmul;
  92.  
  93.  
  94. PROCEDURE rcmul ( a : LONGREAL;  x: LONGCOMPLEX ) : LONGCOMPLEX;
  95. (* returns a*x *)
  96. BEGIN
  97.   INLINE
  98.   (
  99.   PUSH    EDI
  100.   MOV     EDI, DWORD PTR a[ EBP + 8 ]   (* load address of result *)
  101.   FLD     x.r[ EBP ]                    (* r *)
  102.   FLD     x.i[ EBP ]                    (* i r *)
  103.   FLD       a[ EBP ]                    (* a i r *)
  104.   FMUL    ST(2), ST                     (* a i a*r *)
  105.   FMULP   ST(1), ST                     (* a*i a*r *)
  106.   FSTP    QWORD PTR [ EDI + 8 ]         (* rcmul.i := ST *)
  107.   FSTP    QWORD PTR [ EDI ]             (* rcmul.r := ST *)
  108.   FWAIT
  109.   POP     EDI
  110.   );
  111. END rcmul;
  112.  
  113.  
  114. PROCEDURE cdiv( x, y : LONGCOMPLEX ) : LONGCOMPLEX;
  115. (* returns x/y *)
  116. BEGIN
  117.   INLINE
  118.   (
  119.   PUSH    EDI
  120.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  121.   FLD     x.r[ EBP ]                    (* r *)
  122.   FLD     x.i[ EBP ]                    (* i r *)
  123.   FLD     y.r[ EBP ]                    (* rr i r *)
  124.   FLD     y.i[ EBP ]                    (* ii rr i r *)
  125.   FCHS                                  (* -ii rr i r *)
  126.   FLD     ST                            (* -ii -ii rr i r *)
  127.   FMUL    ST, ST(3)                     (* -i*ii -ii rr i r *)
  128.   FXCH    ST(1)                         (* -ii -i*ii rr i r *)
  129.   FMUL    ST, ST(4)                     (* -r*ii -i*ii rr i r *)
  130.   FXCH    ST(2)                         (* rr -i*ii -r*ii i r *)
  131.   FMUL    ST(3), ST                     (* rr -i*ii -r*ii i*rr r *)
  132.   FMULP   ST(4), ST                     (* -i*ii -r*ii i*rr r*rr *)
  133.   FSUBP   ST(3), ST                     (* -r*ii i*rr (r*rr+i*ii) *)
  134.   FADDP   ST(1), ST                     (* (-r*ii+i*rr) (r*rr+i*ii) *)
  135.   FLD     x.r[ EBP ]                    (* x.r y.i y.r *)
  136.   FMUL    ST, ST                        (* x.r*x.r y.i y.r *)
  137.   FLD     x.i[ EBP ]                    (* x.i x.r*x.r y.i y.r *)
  138.   FMUL    ST, ST                        (* x.i*x.i x.r*x.r y.i y.r *)
  139.   FADDP   ST(1), ST                     (* r y.i y.r *)
  140.   FDIV    ST(1), ST                     (* r y.i/r y.r *)
  141.   FDIVP   ST(2), ST                     (* y.i/r y.r/r *)
  142.   FSTP    QWORD PTR [ EDI + 8 ]         (* cmul.i := ST *)
  143.   FSTP    QWORD PTR [ EDI ]             (* cmul.r := ST *)
  144.   FWAIT
  145.   POP     EDI
  146.   );
  147. END cdiv;
  148.  
  149.  
  150. PROCEDURE cdivr ( x: LONGCOMPLEX; a : LONGREAL ) : LONGCOMPLEX;
  151. (* returns x/a *)
  152. BEGIN
  153.   INLINE
  154.   (
  155.   PUSH    EDI
  156.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  157.   FLD     x.r[ EBP ]                    (* r *)
  158.   FLD     x.i[ EBP ]                    (* i r *)
  159.   FLD       a[ EBP ]                    (* a i r *)
  160.   FDIV    ST(2), ST                     (* a i r/a *)
  161.   FDIVP   ST(1), ST                     (* i/a r/a *)
  162.   FSTP    QWORD PTR [ EDI + 8 ]         (* cdivr.i := ST *)
  163.   FSTP    QWORD PTR [ EDI ]             (* cdivr.r := ST *)
  164.   FWAIT
  165.   POP     EDI
  166.   );
  167. END cdivr;
  168.  
  169.  
  170. PROCEDURE csqr( x : LONGCOMPLEX ) : LONGCOMPLEX;
  171. (* returns x*x *)
  172. BEGIN
  173.   INLINE
  174.   (
  175.   PUSH    EDI
  176.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  177.   FLD     x.r[ EBP ]
  178.   FLD     x.i[ EBP ]                    (* i r *)
  179.   FLD     ST                            (* i i r *)
  180.   FMUL    ST(1), ST                     (* i i*i r *)
  181.   FMUL    ST, ST(2)                     (* r*i i*i r *)
  182.   FADD    ST, ST                        (* 2*r*i i*i r *)
  183.   FXCH    ST(2)                         (* r i*i 2*r*i *)
  184.   FMUL    ST, ST                        (* r*r i*i 2*r*i *)
  185.   FSUBRP  ST(1), ST                     (* (r*r - i*i) 2*r*i *)
  186.   FSTP    QWORD PTR [ EDI ]             (* csqr.r := ST *)
  187.   FSTP    QWORD PTR [ EDI + 8 ]         (* csqr.i := ST *)
  188.   FWAIT
  189.   POP     EDI
  190.   );
  191. END csqr;
  192.  
  193.  
  194. PROCEDURE csqrt( x : LONGCOMPLEX ) : LONGCOMPLEX;
  195. (* returns complex square root of x *)
  196. BEGIN
  197.   INLINE
  198.   (
  199.   PUSH    EDI
  200.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* load address of result *)
  201.   FLD     x.r[ EBP ]
  202.   FLD     x.i[ EBP ]                    (* i r *)
  203.   FTST
  204.   FSTSW   AX
  205.   SAHF
  206.   JNE     L1
  207.   (* fall through here if imaginary part is 0 *)
  208.   FXCH    ST(1)                         (* r i *)
  209.   FTST
  210.   FSTSW   AX
  211.   SAHF
  212.   JB      L0
  213.   (* fall through if real part is > 0 *)
  214.   FSQRT                                 (* sqrt(r) i *)
  215.   FXCH    ST(1)                         (* i sqrt(r) *)
  216.   JMP     Done
  217. L0:
  218.   FABS
  219.   FSQRT                                 (* sqrt(abs(r))  i *)
  220.   JMP     Done
  221. L1:   (* come here if imaginary part not zero *)
  222.   FXCH    ST(1)                         (* r i *)
  223.   FTST
  224.   FSTSW   AX
  225.   SAHF
  226.   JNE     L2
  227.   (* fall through here if real part is 0 *)
  228.   FSTP    ST                            (* i *)
  229.   FLD1                                  (* 1 i *)
  230.   FCHS                                  (* -1 i *)
  231.   FXCH    ST(1)                         (* i -1 *)
  232.   FSCALE                                (* i/2 -1 *)
  233.   FSQRT                                 (* sqrt(i/2) -1 *)
  234.   FST     ST(1)                         (* sqrt(i/2) sqrt(i/2) *)
  235.   JMP     Done
  236. L2:   (* come here if both real and imaginary parts not zero *)
  237.   (* 80x87 stack : r i *)
  238.   FLD     ST                            (* r r i *)
  239.   FMUL    ST, ST                        (* r*r r i *)
  240.   FLD     ST(2)                         (* i r*r r i *)
  241.   FMUL    ST, ST                        (* i*i r*r r i *)
  242.   FADDP   ST(1), ST                     (* (r*r+i*i) r i *)
  243.   FSQRT                                 (* sqrt(r*r+i*i) r i *)
  244.   JB      L3
  245.   (* come here if real part > 0 *)
  246.   FADD    ST, ST(1)                     (* r+sqrt(r*r+i*i) r i *)
  247.   FLD1                                  (* 1 r+sqrt(r*r+i*i) r i *)
  248.   FCHS                                  (* -1 r+sqrt(r*r+i*i) r i *)
  249.   FXCH    ST(1)                         (* r+sqrt(r*r+i*i) -1 r i *)
  250.   FSCALE                                (* (r+sqrt(r*r+i*i))/2 -1 r i *)
  251.   FSQRT                                 (* rr -1 r i *)
  252.   FXCH    ST(3)                         (* i -1 r rr *)
  253.   FDIV    ST, ST(3)                     (* i/rr -1 r rr *)
  254.   FSCALE                                (* ii -1 r rr *)
  255.   FSTP    ST(1)                         (* ii r rr *)
  256.   FSTP    ST(1)                         (* ii rr *)
  257.   JMP     Done
  258. L3:   (* come here if real part < 0 *)
  259.   FSUB    ST, ST(1)                     (* -r+sqrt(r*r+i*i) r i *)
  260.   FLD1                                  (* 1 -r+sqrt(r*r+i*i) r i *)
  261.   FCHS                                  (* -1 -r+sqr t(r*r+i*i) r i *)
  262.   FXCH    ST(1)                         (* -r+sqrt(r*r+i*i) -1 r i *)
  263.   FSCALE                                (* (-r+sqrt(r*r+i*i))/2 -1 r i *)
  264.   FSQRT                                 (* ii -1 r i *)
  265.   FXCH    ST(3)                         (* i -1 r ii *)
  266.   FDIV    ST, ST(3)                     (* i/ii -1 r ii *)
  267.   FSCALE                                (* rr -1 r ii *)
  268.   FSTP    ST(1)                         (* rr r ii *)
  269.   FSTP    ST(1)                         (* rr ii *)
  270.   FXCH    ST(1)                         (* ii rr *)
  271. Done:
  272.   FSTP    QWORD PTR [ EDI + 8 ]         (* csqrt.i := ST *)
  273.   FSTP    QWORD PTR [ EDI ]             (* csqrt.r := ST *)
  274.   FWAIT
  275.   POP     EDI
  276.   );
  277. END csqrt;
  278.  
  279.  
  280. PROCEDURE lcmplx( r, i : LONGREAL ) : LONGCOMPLEX;
  281. (* returns complex number *)
  282. BEGIN
  283.   INLINE
  284.   (
  285.   PUSH    EDI
  286.   MOV     EDI, DWORD PTR r[ EBP + 8 ]  (* hidden address of result *)
  287.   FLD     r[ EBP ]                     (* real part *)
  288.   FSTP    QWORD PTR [ EDI ]            (* lcmplx.r := ST *)
  289.   FLD     i[ EBP ]                     (* imaginary part *)
  290.   FSTP    QWORD PTR [ EDI + 8 ]        (* lcmplx.i := ST *)
  291.   FWAIT
  292.   POP     EDI
  293.   );
  294. END lcmplx;
  295.  
  296.  
  297. PROCEDURE cmag( x : LONGCOMPLEX ) : LONGREAL;
  298. (* returns magnitude of x = sqrt(x.r*x.r+x.i*x.i) *)
  299. BEGIN
  300.   INLINE
  301.   (
  302.   FLD     x.r[ EBP ]
  303.   FMUL    ST, ST
  304.   FSTP    x.r[ EBP ]
  305.   FLD     x.i[ EBP ]
  306.   FMUL    ST, ST
  307.   FADD    x.r[ EBP ]
  308.   FSQRT
  309.   );
  310. END cmag;
  311.  
  312.  
  313. PROCEDURE conjg ( x : LONGCOMPLEX ) : LONGCOMPLEX;
  314. (* returns complex conjugate of x *)
  315. BEGIN
  316.   INLINE
  317.   (
  318.   PUSH    EDI
  319.   MOV     EDI, DWORD PTR x[ EBP + 16 ]  (* hidden address of result *)
  320.   FLD     x.r[ EBP ]                    (* real part *)
  321.   FSTP    QWORD PTR [ EDI ]             (* conjg.r := ST *)
  322.   FLD     x.i[ EBP ]                    (* imaginary part *)
  323.   FCHS
  324.   FSTP    QWORD PTR [ EDI + 8 ]         (* conjg.i := ST *)
  325.   FWAIT
  326.   POP     EDI
  327.   );
  328. END conjg;
  329.  
  330.  
  331. PROCEDURE cexp ( x : LONGCOMPLEX ) : LONGCOMPLEX;
  332. (* returns complex exp(x) *)
  333. VAR
  334.   ex     : LONGREAL;
  335. BEGIN
  336.   ex := MathLib0.exp( x.r );
  337.   INLINE
  338.   (
  339.   PUSH    EDI
  340.   MOV     EDI, DWORD PTR x[ EBP + 16 ] (* hidden address of result *)
  341.   FLD     x.i[ EBP ]                   (* x.i *)
  342.   FSINCOS                              (* cos(x.i) sin(x.i) *)
  343.   FMUL    ex[ EBP ]                    (* ex*cos(x.i) sin(x.i) *)
  344.   FSTP    QWORD PTR [ EDI ]            (* cexp.r := ST *)
  345.   FMUL    ex[ EBP ]                    (* ex*sin(x.i) *)
  346.   FSTP    QWORD PTR [ EDI + 8 ]        (* cexp.i := ST *)
  347.   FWAIT
  348.   POP     EDI
  349.   );
  350. END cexp;
  351.  
  352.  
  353. PROCEDURE cln ( x : LONGCOMPLEX ) : LONGCOMPLEX;
  354. (* returns complex ln(x) *)
  355. BEGIN
  356.   INLINE
  357.   (
  358.   PUSH    EDI
  359.   MOV     EDI, DWORD PTR x[ EBP + 16 ] (* hidden address of result *)
  360.   FLD     x.r[ EBP ]                   (* r *)
  361.   FLD     x.i[ EBP ]                   (* i r *)
  362.   FLD1                                 (* 1 i r *)
  363.   FCHS                                 (* -1 i r *)
  364.   FLD     ST(1)                        (* i -1 i r *)
  365.   FMUL    ST, ST                       (* i*i -1 i r *)
  366.   FLD     ST(3)                        (* r i*i -1 i r *)
  367.   FMUL    ST, ST                       (* r*r i*i -1 i r *)
  368.   FADDP   ST(1), ST                    (* (r*r+i*i) -1 i r *)
  369.   FLDLN2                               (* log2e (r*r+i*i) -1 i r *)
  370.   FXCH    ST(1)                        (* (r*r+i*i)/2 log2e -1 i r *)
  371.   FYL2X                                (* ln(r*r+i*i) -1 i r *)
  372.   FSCALE
  373.   FSTP    ST(1)                        (* 0.5*ln(r*r+i*i) i r *)
  374.   FSTP    QWORD PTR [ EDI ]            (* cln.r := ST *)
  375.   FXCH    ST(1)                        (* r i *)
  376.   FPATAN                               (* atan(ST(1)/ST) *)
  377.   FSTP    QWORD PTR [ EDI + 8 ]        (* cln.i := ST *)
  378.   FWAIT
  379.   POP     EDI
  380.   );
  381. END cln;
  382.  
  383.  
  384. PROCEDURE csin ( x : LONGCOMPLEX ) : LONGCOMPLEX;
  385. (* returns complex sin(x) *)
  386.   ex     : LONGREAL;
  387. BEGIN
  388.   ex := MathLib0.exp( x.i );
  389.   INLINE
  390.   (
  391.   PUSH    EDI
  392.   MOV     EDI, DWORD PTR x[ EBP + 16 ] (* hidden address of result *)
  393.   FLD     x.r[ EBP ]                   (* x.r *)
  394.   FSINCOS                              (* cos(x.r) sin(x.r) *)
  395.   FLD1                                 (* 1 c s *)
  396.   FLD     ex[ EBP ]                    (* ex 1 c s *)
  397.   FDIV    ST(1), ST                    (* ex 1/ex c s *)
  398.   FLD     ST                           (* ex ex 1/ex c s *)
  399.   FADD    ST, ST(2)                    (* (ex+1/ex) ex 1/ex c s *)
  400.   FMULP   ST(4), ST                    (* ex 1/ex c (ex+1/ex)*s *)
  401.   FSUBRP  ST(1), ST                    (* (ex-1/ex) c (ex+1/ex)*s *)
  402.   FMULP   ST(1), ST                    (* (ex-1/ex)*c (ex+1/ex)*s *)
  403.   FLD1
  404.   FCHS                                 (* -1 (ex-1/ex)*c (ex+1/ex)*s *)
  405.   FXCH    ST(1)                        (* (ex-1/ex)*c -1 (ex+1/ex)*s *)
  406.   FSCALE
  407.   FSTP    QWORD PTR [ EDI + 8 ]        (* csin.i := ST *)
  408.   FXCH    ST(1)                        (* (ex+1/ex)*s -1 *)
  409.   FSCALE
  410.   FSTP    ST(1)
  411.   FSTP    QWORD PTR [ EDI ]
  412.   FWAIT
  413.   POP     EDI
  414.   );
  415. END csin;
  416.  
  417.  
  418. PROCEDURE ccos ( x : LONGCOMPLEX ) : LONGCOMPLEX;
  419. (* returns complex cos(x) *)
  420.   ex     : LONGREAL;
  421. BEGIN
  422.   ex := MathLib0.exp( x.i );
  423.   INLINE
  424.   (
  425.   PUSH    EDI
  426.   MOV     EDI, DWORD PTR x[ EBP + 16 ] (* hidden address of result *)
  427.   FLD     x.r[ EBP ]                   (* x.r *)
  428.   FSINCOS                              (* cos(x.r) sin(x.r) *)
  429.   FLD1                                 (* 1 c s *)
  430.   FLD     ex[ EBP ]                    (* ex 1 c s *)
  431.   FDIV    ST(1), ST                    (* ex 1/ex c s *)
  432.   FLD     ST                           (* ex ex 1/ex c s *)
  433.   FADD    ST, ST(2)                    (* (ex+1/ex) ex 1/ex c s *)
  434.   FMULP   ST(3), ST                    (* ex 1/ex (ex+1/ex)*c s *)
  435.   FSUBP   ST(1), ST                    (* (1/ex-ex) (ex+1/ex)*c s *)
  436.   FMULP   ST(2), ST                    (* (ex+1/ex)*c (1/ex-ex)*s *)
  437.   FLD1
  438.   FCHS
  439.   FXCH    ST(1)                        (* (ex+1/ex)*c -1 (1/ex-ex)*s *)
  440.   FSCALE
  441.   FSTP    QWORD PTR [ EDI ]            (* ccos.r := ST *)
  442.   FXCH    ST(1)                        (* (1/ex-ex)*s -1 *)
  443.   FSCALE
  444.   FSTP    QWORD PTR [ EDI + 8 ]        (* ccos.i := ST *)
  445.   FWAIT
  446.   POP     EDI
  447.   );
  448. END ccos;
  449.  
  450.  
  451. PROCEDURE WriteLongComplex( x : LONGCOMPLEX; n : INTEGER );
  452. (* write complex number x *)
  453. BEGIN
  454.   InOut.Write('(');
  455.   RealInOut.WriteLongReal( x.r, n );
  456.   InOut.Write(',');
  457.   RealInOut.WriteLongReal( x.i, n );
  458.   InOut.Write(')');
  459. END WriteLongComplex;
  460.  
  461.  
  462. BEGIN
  463. END ComplexLib.
  464.  
  465.