home *** CD-ROM | disk | FTP | other *** search
/ Die Ultimative Software-P…i Collection 1996 & 1997 / Die Ultimative Software-Pakete CD-ROM fur Atari Collection 1996 & 1997.iso / g / gnu_c / pmlsrc23.zoo / pmlsrc / cdiv.c < prev    next >
Encoding:
C/C++ Source or Header  |  1994-03-19  |  5.6 KB  |  269 lines
  1. /************************************************************************
  2.  *                                    *
  3.  *                N O T I C E                *
  4.  *                                    *
  5.  *            Copyright Abandoned, 1987, Fred Fish        *
  6.  *                                    *
  7.  *    This previously copyrighted work has been placed into the    *
  8.  *    public domain by the author (Fred Fish) and may be freely used    *
  9.  *    for any purpose, private or commercial.  I would appreciate    *
  10.  *    it, as a courtesy, if this notice is left in all copies and    *
  11.  *    derivative works.  Thank you, and enjoy...            *
  12.  *                                    *
  13.  *    The author makes no warranty of any kind with respect to this    *
  14.  *    product and explicitly disclaims any implied warranties of    *
  15.  *    merchantability or fitness for any particular purpose.        *
  16.  *                                    *
  17.  ************************************************************************
  18.  */
  19.  
  20.  
  21. /*
  22.  *  FUNCTION
  23.  *
  24.  *    cdiv   double precision complex division
  25.  *
  26.  *  KEY WORDS
  27.  *
  28.  *    cdiv
  29.  *    complex functions
  30.  *    machine independent routines
  31.  *    math libraries
  32.  *
  33.  *  DESCRIPTION
  34.  *
  35.  *    Computes double precision complex result of division of
  36.  *    first double precision complex argument by second double
  37.  *    precision complex argument.
  38.  *
  39.  *  USAGE
  40.  *
  41.  *    COMPLEX cdiv (numerator, denominator)
  42.  *    COMPLEX numerator;
  43.  *    COMPLEX denominator;
  44.  *
  45.  *  PROGRAMMER
  46.  *
  47.  *    Fred Fish
  48.  *    Tempe, Az 85281
  49.  *    (602) 966-8871
  50.  *
  51.  *  INTERNALS
  52.  *
  53.  *    Computes cdiv(znum,zden) from:
  54.  *
  55.  *        1.    Let znum = a + j b
  56.  *            Let zden = c + j d
  57.  *
  58.  *        2.    denom = c*c + d*d
  59.  *
  60.  *        3.    If denom is zero then log error,
  61.  *            set r_cdiv = maximum floating value,
  62.  *            i_cdiv = 0, and go to step 5.
  63.  *
  64.  *        4.    r_cdiv = (a*c + b*d) / denom
  65.  *            i_cdiv = (c*b - a*d) / denom
  66.  *
  67.  *        5.    Then cdiv(znum,zden) = r_cdiv + j i_cdiv
  68.  *
  69.  */
  70.  
  71. #if !defined (__M68881__) && !defined (sfp004)
  72.  
  73. #include <stdio.h>
  74. #include <math.h>
  75. #include "pml.h"
  76.  
  77. COMPLEX cdiv (znum, zden)
  78. COMPLEX znum;
  79. COMPLEX zden;
  80. {
  81.     COMPLEX result;
  82.     double denom;
  83.     struct exception xcpt;
  84.  
  85.     result.real = ((znum.real*zden.real)+(znum.imag*zden.imag));
  86.     result.imag = ((zden.real*znum.imag)-(znum.real*zden.imag));
  87.     denom = (zden.real * zden.real) + (zden.imag * zden.imag);
  88.  
  89.     if (denom == 0.0) {
  90. #ifdef    ERROR_CHECK
  91.     xcpt.type = SING;
  92.     xcpt.name = "cdiv";
  93.     xcpt.arg1 = denom;
  94.     if (!matherr (&xcpt)) {
  95.         fprintf (stderr, "%s:  ZERO_CMPLX_DENOMINATOR \n", xcpt.name);
  96.         xcpt.retval = 0.0;    /* useless in this context */
  97.         errno = ERANGE;    /* should be EDOM if real or imag == 0    */
  98.     }
  99.     if( result.real >= 0.0)    result.real = HUGE_VAL;    
  100.                     /* still wrong, == 0 should yield NAN */
  101.     else            result.real = -HUGE_VAL;    
  102.     if( result.imag >= 0.0) result.imag = HUGE_VAL;
  103.                     /* still wrong, == 0 should yield NAN */
  104.     else            result.imag = -HUGE_VAL;    
  105. #else    ERROR_CHECK
  106.     result.real /= denom;
  107.     result.imag /= denom;
  108. #endif    ERROR_CHECK
  109.     } else {
  110.     result.real /= denom;
  111.     result.imag /= denom;
  112.     }
  113.     return (result);
  114. }
  115. #endif !defined (__M68881__) && !defined (sfp004)
  116. #ifdef    __M68881__
  117. __asm("
  118. .text
  119. .even
  120. _funcname:
  121.     .ascii    \"cdiv\\0\"
  122.     .even
  123.  
  124. .globl    _cdiv
  125. _cdiv:
  126.     fmoved    sp@(4),fp0
  127.     fmoved    sp@(12),fp1
  128.     fmoved    sp@(20),fp2
  129.     fmoved    sp@(28),fp3
  130.     fmovex    fp0,fp4
  131.     movel    a1,d0        | pointer to result
  132.  
  133.     fmovex    fp2,fp5
  134.     fmulx    fp2,fp5
  135.     fmovex    fp3,fp6
  136.     fmulx    fp3,fp6
  137.     faddx    fp6,fp5
  138.  
  139.     fmulx    fp2,fp4
  140.     fmulx    fp3,fp0
  141.     fmulx    fp1,fp2    
  142.     fmulx    fp1,fp3
  143.     faddx    fp3,fp4
  144.     fdivx    fp5,fp4
  145.     fsubx    fp0,fp2
  146.     fdivx    fp5,fp2
  147.  
  148.     fmoved    fp4,a1@
  149.     fmoved    fp2,a1@(8)
  150. ");    /* end asm    */
  151. #endif    __M68881__
  152.  
  153. #ifdef    sfp004
  154. __asm("
  155.  
  156. comm =     -6
  157. resp =    -16
  158. zahl =      0
  159.  
  160. .even
  161. .text
  162. .even
  163. _funcname:
  164.     .ascii    \"cdiv\\0\"
  165.     .even
  166. .text
  167. .even
  168. .globl    _cdiv
  169. _cdiv:
  170.  
  171.     lea    0xfffa50,a0
  172.  
  173.     movew    #0x5400,a0@(comm)    | z1.real -> fp0
  174.     movel    a1,d0            | pointer to result
  175.     .long    0x0c688900, 0xfff067f8
  176.     movel    a7@(4),a0@        | load arg_hi
  177.     movel    a7@(8),a0@        | load arg_low
  178.  
  179.     movew    #0x5480,a0@(comm)    | z1.imag -> fp1
  180.     .long    0x0c688900, 0xfff067f8
  181.     movel    a7@(12),a0@        | load arg_hi
  182.     movel    a7@(16),a0@        | load arg_low
  183.  
  184.     movew    #0x5500,a0@(comm)    | z2.real -> fp2
  185.     .long    0x0c688900, 0xfff067f8
  186.     movel    a7@(20),a0@        | load arg_hi
  187.     movel    a7@(24),a0@        | load arg_low
  188.  
  189.     movew    #0x5580,a0@(comm)    | z2.imag -> fp3
  190.     .long    0x0c688900, 0xfff067f8
  191.     movel    a7@(28),a0@        | load arg_hi
  192.     movel    a7@(32),a0@        | load arg_low
  193.  
  194.     movew    #0x0200,a0@(comm)    | copy fp0 to fp4
  195.     .word    0x4a68,0xfff0,0x6bfa    | test
  196.  
  197.  
  198. |    fmovex    fp2,fp5
  199.     movew    #0x0a80,a0@(comm)
  200.     .word    0x4a68,0xfff0,0x6bfa    | test
  201. |    fmulx    fp2,fp5
  202.     movew    #0x0aa3,a0@(comm)
  203.     .word    0x4a68,0xfff0,0x6bfa    | test
  204. |    fmovex    fp3,fp6
  205.     movew    #0x0f00,a0@(comm)
  206.     .word    0x4a68,0xfff0,0x6bfa    | test
  207. |    fmulx    fp3,fp6
  208.     movew    #0x0f23,a0@(comm)
  209.     .word    0x4a68,0xfff0,0x6bfa    | test
  210. |    faddx    fp6,fp5
  211.     movew    #0x1aa2,a0@(comm)
  212.     .word    0x4a68,0xfff0,0x6bfa    | test
  213.  
  214. |    fmulx    fp2,fp4
  215.     movew    #0x0a23,a0@(comm)
  216.     .word    0x4a68,0xfff0,0x6bfa    | test
  217. |    fmulx    fp3,fp0
  218.     movew    #0x0c23,a0@(comm)
  219.     .word    0x4a68,0xfff0,0x6bfa    | test
  220. |    fmulx    fp1,fp2
  221.     movew    #0x0523,a0@(comm)
  222.     .word    0x4a68,0xfff0,0x6bfa    | test
  223. |    fmulx    fp1,fp3
  224.     movew    #0x05a3,a0@(comm)
  225.     .word    0x4a68,0xfff0,0x6bfa    | test
  226.  
  227. |    faddx    fp3,fp4
  228.     movew    #0x0e22,a0@(comm)
  229.     .word    0x4a68,0xfff0,0x6bfa    | test
  230. |    fdivx    fp5,fp4
  231.     movew    #0x1620,a0@(comm)
  232.     .word    0x4a68,0xfff0,0x6bfa    | test
  233. |    fsubx    fp0,fp2
  234.     movew    #0x0128,a0@(comm)
  235.     .word    0x4a68,0xfff0,0x6bfa    | test
  236. |    fdivx    fp5,fp2
  237.     movew    #0x1520,a0@(comm)
  238.     .word    0x4a68,0xfff0,0x6bfa    | test
  239.  
  240.  
  241. |    fmoved    fp4,a1@
  242.     movew    #0x7600,a0@(comm)        | 
  243.     .long    0x0c688900, 0xfff067f8
  244.     movel    a0@,a1@
  245.     movel    a0@,a1@(4)
  246.  
  247. |    fmoved    fp2,a1@(8)
  248.     movew    #0x7500,a0@(comm)        | 
  249.     .long    0x0c688900, 0xfff067f8
  250.     movel    a0@,a1@(8)
  251.     movel    a0@,a1@(12)
  252. ");    /* end asm    */
  253. #endif    sfp004
  254.  
  255. #if defined (__M68881__) || defined (sfp004)
  256. # ifdef ERROR_CHECK    /* no error checking for now    */
  257.  
  258. __asm("    
  259.     pea    _funcname
  260.     jmp    c_err_check
  261. ");    /* end asm    */
  262.  
  263. # else  ERROR_CHECK
  264.  
  265. __asm("rts");
  266.  
  267. # endif ERROR_CHECK
  268. #endif defined (__M68881__) || defined (sfp004)
  269.