home *** CD-ROM | disk | FTP | other *** search
/ Club Amiga de Montreal - CAM / CAM_CD_1.iso / files / 221.lha / src / complex.d < prev    next >
Text File  |  1996-02-15  |  9KB  |  405 lines
  1. #drinc:float.g
  2. #drinc:util.g
  3.  
  4. type Complex_t = struct {
  5.     float cmpl_real;
  6.     float cmpl_imag;
  7. };
  8.  
  9. uint STACK_HEIGHT = 20;
  10.  
  11. [STACK_HEIGHT] Complex_t Stack;
  12.  
  13. *Complex_t StackTop;
  14.  
  15. /*
  16.  * cmplInitialize - set up the world.
  17.  */
  18.  
  19. proc cmplInitialize()void:
  20.  
  21.     StackTop := &Stack[STACK_HEIGHT];
  22. corp;
  23.  
  24. /*
  25.  * _cmplpsh - calls generated by the compiler to push complex values onto our
  26.  *    internal stack. The address of the complex value will be on the normal
  27.  *    stack as a parameter.
  28.  */
  29.  
  30. proc _cmplpsh(*Complex_t n)void:
  31.  
  32.     StackTop := StackTop - sizeof(Complex_t);
  33.     StackTop* := n*;
  34. corp;
  35.  
  36. /*
  37.  * _cmplpop - calls generated by the compiler to pop complex values from our
  38.  *    internal stack into a destination variable. The address of the
  39.  *    destination variable will be on the normal stack as a parameter.
  40.  */
  41.  
  42. proc _cmplpop(*Complex_t n)void:
  43.  
  44.     n* := StackTop*;
  45.     StackTop := StackTop + sizeof(Complex_t);
  46. corp;
  47.  
  48. /*
  49.  * _cmplput - calls generated by the compiler to write a complex value in
  50.  *    text mode.
  51.  */
  52.  
  53. proc _cmplput()void:
  54.  
  55.     write(*; '(', StackTop*.cmpl_real, ',', StackTop*.cmpl_imag, ')');
  56.     StackTop := StackTop + sizeof(Complex_t);
  57. corp;
  58.  
  59. /*
  60.  * _cmplget - calls generated by the compiler to read a complex value in
  61.  *    text mode. Various error conditions can occur. An option here would
  62.  *    be to make the parentheses and the comma required. The value read in
  63.  *    is left on our internal stack.
  64.  */
  65.  
  66. proc _cmplget()void:
  67.     extern
  68.     _d_channelSkip()void,
  69.     _d_channelGetChar()char,
  70.     _d_channelUnGetChar(char ch)void,
  71.     _d_channelHadError()bool,
  72.     _d_channelError(ushort err)void;
  73.     register *Complex_t p;
  74.     register char ch;
  75.     bool gotParen;
  76.  
  77.     p := StackTop - sizeof(Complex_t);
  78.     StackTop := p;
  79.     _d_channelSkip();
  80.     ch := _d_channelGetChar();
  81.     if ch = '(' then
  82.     gotParen := true;
  83.     else
  84.     _d_channelUnGetChar(ch);
  85.     gotParen := false;
  86.     fi;
  87.     read(*; p*.cmpl_real);
  88.     if not _d_channelHadError() then
  89.     _d_channelSkip();
  90.     ch := _d_channelGetChar();
  91.     if ch ~= ',' then
  92.         _d_channelUnGetChar(ch);
  93.     fi;
  94.     read(*; p*.cmpl_imag);
  95.     if not _d_channelHadError() then
  96.         if gotParen then
  97.         _d_channelSkip();
  98.         ch := _d_channelGetChar();
  99.         if ch ~= ')' then
  100.             _d_channelUnGetChar(ch);
  101.             _d_channelError(CH_BADCHAR);
  102.         fi;
  103.         fi;
  104.     fi;
  105.     fi;
  106. corp;
  107.  
  108. /*
  109.  * _cmplneg - compiler generated for the '-' unary operator.
  110.  */
  111.  
  112. proc _cmplneg()void:
  113.     register *Complex_t p;
  114.  
  115.     p := StackTop;
  116.     p*.cmpl_real := - p*.cmpl_real;
  117.     p*.cmpl_imag := - p*.cmpl_imag;
  118. corp;
  119.  
  120. /*
  121.  * cmpladd - compiler generated for the '+' binary operator.
  122.  */
  123.  
  124. proc _cmpladd()void:
  125.     register *Complex_t stackTop;
  126.  
  127.     stackTop := StackTop;
  128.     (stackTop + sizeof(Complex_t))*.cmpl_real :=
  129.     (stackTop + sizeof(Complex_t))*.cmpl_real + stackTop*.cmpl_real;
  130.     (stackTop + sizeof(Complex_t))*.cmpl_imag :=
  131.     (stackTop + sizeof(Complex_t))*.cmpl_imag + stackTop*.cmpl_imag;
  132.     StackTop := stackTop + sizeof(Complex_t);
  133. corp;
  134.  
  135. /*
  136.  * _cmplsub - compler generated for the '-' binary operator.
  137.  */
  138.  
  139. proc _cmplsub()void:
  140.     register *Complex_t stackTop;
  141.  
  142.     stackTop := StackTop;
  143.     (stackTop + sizeof(Complex_t))*.cmpl_real :=
  144.     (stackTop + sizeof(Complex_t))*.cmpl_real - stackTop*.cmpl_real;
  145.     (stackTop + sizeof(Complex_t))*.cmpl_imag :=
  146.     (stackTop + sizeof(Complex_t))*.cmpl_imag - stackTop*.cmpl_imag;
  147.     StackTop := stackTop + sizeof(Complex_t);
  148. corp;
  149.  
  150. /*
  151.  * _cmplmul - compiler generated for the '*' binary operator.
  152.  */
  153.  
  154. proc _cmplmul()void:
  155.     register *Complex_t l, r;
  156.     float temp;
  157.  
  158.     r := StackTop;
  159.     l := r + sizeof(Complex_t);
  160.     temp := l*.cmpl_real;
  161.     l*.cmpl_real := temp * r*.cmpl_real - l*.cmpl_imag * r*.cmpl_imag;
  162.     l*.cmpl_imag := l*.cmpl_imag * r*.cmpl_real + temp * r*.cmpl_imag;
  163.     StackTop := l;
  164. corp;
  165.  
  166. /*
  167.  * _cmpldiv - compiler generated for the '/' binary operator.
  168.  */
  169.  
  170. proc _cmpldiv()void:
  171.     register *Complex_t l, r;
  172.     float t1, t2;
  173.  
  174.     r := StackTop;
  175.     l := r + sizeof(Complex_t);
  176.     t1 := l*.cmpl_real;
  177.     t2 := 1.0 / (r*.cmpl_real * r*.cmpl_real + r*.cmpl_imag * r*.cmpl_imag);
  178.     l*.cmpl_real := (t1 * r*.cmpl_real + l*.cmpl_imag * r*.cmpl_imag) * t2;
  179.     l*.cmpl_imag := (l*.cmpl_imag * r*.cmpl_real - t1 * r*.cmpl_imag) * t2;
  180.     StackTop := l;
  181. corp;
  182.  
  183. /*****************************************************************************\
  184. *                                          *
  185. * The routines from here on down could have been in another source file, in   *
  186. * which case they would have had one or more 'complex' arguments or result,   *
  187. * and would not explicitly reference 'StackTop'. This would have made them    *
  188. * slower, however, since they would have been calling '_cmplpsh', '_cmplpop', *
  189. * etc., so they are included here to allow them direct access. All of these   *
  190. * routines are called directly by the user, and not by compiler-generated     *
  191. * code.                                       *
  192. *                                          *
  193. \*****************************************************************************/
  194.  
  195. /*
  196.  * cmpl - return a complex build from two float values.
  197.  *    This differs from the compiler's ability to make a complex constant
  198.  *    from two float constants in that the expressions need not be
  199.  *    constants, and the construction is done at run-time.
  200.  */
  201.  
  202. proc cmpl(float re, im)void:
  203.  
  204.     StackTop := StackTop - sizeof(Complex_t);
  205.     StackTop*.cmpl_real := re;
  206.     StackTop*.cmpl_imag := im;
  207. corp;
  208.  
  209. /*
  210.  * cmplRe - return the real part of a complex number.
  211.  */
  212.  
  213. proc cmplRe()float:
  214.  
  215.     StackTop := StackTop + sizeof(Complex_t);
  216.     (StackTop - sizeof(Complex_t))*.cmpl_real
  217. corp;
  218.  
  219. /*
  220.  * cmplIm - return the imaginary part of a complex number.
  221.  */
  222.  
  223. proc cmplIm()float:
  224.  
  225.     StackTop := StackTop + sizeof(Complex_t);
  226.     (StackTop - sizeof(Complex_t))*.cmpl_imag
  227. corp;
  228.  
  229. /*
  230.  * cmplConj - return the complex conjugate of a value.
  231.  */
  232.  
  233. proc cmplConj()void:
  234.  
  235.     StackTop*.cmpl_imag := - StackTop*.cmpl_imag;
  236. corp;
  237.  
  238. /*
  239.  * cmplNorm - return the complex norm (a float value) of a complex value.
  240.  */
  241.  
  242. proc cmplNorm()float:
  243.     register *Complex_t p;
  244.  
  245.     p := StackTop;
  246.     StackTop := p + sizeof(Complex_t);
  247.     p*.cmpl_real * p*.cmpl_real + p*.cmpl_imag * p*.cmpl_imag
  248. corp;
  249.  
  250. /*
  251.  * cmplAbs - complex absolute value (a float value) of a complex value.
  252.  */
  253.  
  254. proc cmplAbs()float:
  255.  
  256.     sqrt(cmplNorm())
  257. corp;
  258.  
  259. /*
  260.  * cmplSqrt - complex square root of a complex value.
  261.  */
  262.  
  263. proc cmplSqrt()void:
  264.     register *Complex_t p;
  265.     float temp;
  266.  
  267.     p := StackTop;
  268.     temp := cmplAbs();
  269.     StackTop := p;
  270.     p*.cmpl_imag := sqrt((temp - p*.cmpl_real) * 0.5);
  271.     p*.cmpl_real := sqrt((temp + p*.cmpl_real) * 0.5);
  272. corp;
  273.  
  274. /*
  275.  * cmplArg
  276.  */
  277.  
  278. proc cmplArg()float:
  279.     register *Complex_t p;
  280.     float temp;
  281.  
  282.     p := StackTop;
  283.     StackTop := p + sizeof(Complex_t);
  284.     if p*.cmpl_real = 0.0 then
  285.     if p*.cmpl_imag > 0.0 then
  286.         HALF_PI
  287.     elif p*.cmpl_imag < 0.0 then
  288.         - HALF_PI
  289.     else
  290.         0.0
  291.     fi
  292.     else
  293.     temp := atan(p*.cmpl_imag / p*.cmpl_real);
  294.     if p*.cmpl_real > 0.0 then
  295.         temp
  296.     else
  297.         temp + PI
  298.     fi
  299.     fi
  300. corp;
  301.  
  302. /*
  303.  * cmplExp
  304.  */
  305.  
  306. proc cmplExp()void:
  307.     register *Complex_t p;
  308.     float temp;
  309.  
  310.     p := StackTop;
  311.     temp := exp(p*.cmpl_real);
  312.     p*.cmpl_real := temp * cos(p*.cmpl_imag);
  313.     p*.cmpl_imag := temp * sin(p*.cmpl_imag);
  314. corp;
  315.  
  316. /*
  317.  * cmplLog
  318.  */
  319.  
  320. proc cmplLog()void:
  321.     register *Complex_t p;
  322.     float temp;
  323.  
  324.     p := StackTop;
  325.     temp := log(cmplNorm()) * 0.5;
  326.     StackTop := p;
  327.     p*.cmpl_imag := cmplArg();
  328.     StackTop := p;
  329.     p*.cmpl_real := temp;
  330. corp;
  331.  
  332. /*
  333.  * cmplSin
  334.  */
  335.  
  336. proc cmplSin()void:
  337.     register *Complex_t p;
  338.     float temp;
  339.  
  340.     p := StackTop;
  341.     temp := p*.cmpl_real;
  342.     p*.cmpl_real := sin(temp) * cosh(p*.cmpl_imag);
  343.     p*.cmpl_imag := cos(temp) * sinh(p*.cmpl_imag);
  344. corp;
  345.  
  346. /*
  347.  * cmplCos
  348.  */
  349.  
  350. proc cmplCos()void:
  351.     register *Complex_t p;
  352.     float temp;
  353.  
  354.     p := StackTop;
  355.     temp := p*.cmpl_real;
  356.     p*.cmpl_real := cos(temp) * cosh(p*.cmpl_imag);
  357.     p*.cmpl_imag := - (sin(temp) * sinh(p*.cmpl_imag));
  358. corp;
  359.  
  360. /*
  361. proc cmplTan()void:corp;
  362. proc cmplAsin()void:corp;
  363. proc cmplAcos()void:corp;
  364. */
  365.  
  366. /*
  367.  * atanh - utility routine needed below.
  368.  */
  369.  
  370. proc atanh(float n)float:
  371.  
  372.     0.5 * log((1.0 + n) / (1.0 - n))
  373. corp;
  374.  
  375. /*
  376.  * cmplAtan
  377.  */
  378.  
  379. proc cmplAtan()void:
  380.     register *Complex_t p;
  381.     float anorm, snorm;
  382.  
  383.     p := StackTop;
  384.     if p*.cmpl_imag = 0.0 then
  385.     p*.cmpl_real := atan(p*.cmpl_real);
  386.     elif p*.cmpl_real = 0.0 then
  387.     if p*.cmpl_imag <= -1.0 then
  388.         p*.cmpl_imag := atanh(1.0 / p*.cmpl_imag);
  389.         p*.cmpl_real := - HALF_PI;
  390.     elif p*.cmpl_imag >= 1.0 then
  391.         p*.cmpl_imag := atanh(1.0 / p*.cmpl_imag);
  392.         p*.cmpl_real := HALF_PI;
  393.     else
  394.         p*.cmpl_imag := atanh(p*.cmpl_imag);
  395.     fi;
  396.     else
  397.     snorm := cmplNorm();
  398.     StackTop := p;
  399.     anorm := snorm - 1.0;
  400.     anorm := sqrt(anorm * anorm + 4.0 * p*.cmpl_real * p*.cmpl_real);
  401.     p*.cmpl_real := atan(0.5 * (snorm + anorm - 1.0) / p*.cmpl_real);
  402.     p*.cmpl_imag := atanh(0.5 * (snorm - anorm + 1.0) / p*.cmpl_imag);
  403.     fi;
  404. corp;
  405.