home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 3 / 3137 / dflonum.c next >
Encoding:
C/C++ Source or Header  |  1991-03-27  |  22.1 KB  |  1,046 lines
  1. /* A working Gnulib68k, thanks to Scott McCauley for the origonal
  2.    version, and John Dunning (jrd@STONY-BROOK.SCRC.Symbolics.COM)
  3.    who got this working.
  4.    
  5.    Please not that only the double float. stuff is picked up from
  6.    here, the other long-int and single float stuff come from
  7.    fixnum.s and sflonum.s (see flonum.h for the appr. #defs).
  8.    ++jrb
  9.    */
  10.  
  11. /* Subroutines needed by GCC output code for the 68000/20. 
  12.    
  13.    Compile using -O flag with gcc. 
  14.    Use the -m68000 flag if you have a 68000
  15.    
  16.    This package includes 32 bit integer and 64-bit floating
  17.    point routines.
  18.    
  19.    WARNING: alpha-test version (July 1988) -- probably full of bugs.
  20.    If you find a bug, send bugs reports to jsm@phoenix.princeton.edu,
  21.    or
  22.    
  23.    Scott McCauley
  24.    PPPL P. O. Box 451
  25.    Princeton NJ 08543
  26.    
  27.    Known bugs/features:
  28.    
  29.    1) Depending on the version of GCC, this may produce a lot of
  30.    warnings about conversion clashes, etc. It should still compile
  31.    as is. Also, there appears to be a bug in the register-allocation
  32.    parts of gcc that makes the compiler sometimes core dump
  33.    or run into an infinite loop. This version works -- if you
  34.    decide to change routines these compiler bugs can bite very hard....
  35.    
  36.    2) all single-precision gets converted to double precision in the
  37.    math routines (in accordance with K&R), so doing things in
  38.    double precision is faster than single precision....
  39.  
  40. (not any more -- jrd )
  41.    
  42.    3) double precision floating point division may not be accurate to
  43.    the last four or so bits. Most other routines round to the
  44.    lsb.
  45.  
  46. (not any more -- jrd )
  47.    
  48.    4) no support of NaN and Inf
  49.    
  50.    5) beware of undefined operations: i.e. a float to integer conversion
  51.    when the float is larger than MAXIINT.
  52.    
  53.    */
  54.  
  55. #include "flonum.h"
  56.  
  57. #ifdef __GCC_OLD__
  58. #define umultl _umulsi
  59. #define multl _mulsi3
  60. #define udivl _udivsi3 
  61. #define divl _divsi3
  62. #define ddiv _divdf3
  63. #define qmult _qmult
  64. #define dmult _muldf3
  65. #define dneg _negdf2
  66. #define dadd _adddf3
  67. #define dcmp _cmpdf2
  68. #define dtoul _fixunsdfsi
  69. #define dtol _fixdfsi
  70. #define ltod _floatsidf
  71. #else
  72. #define umultl __umulsi
  73. #define multl __mulsi3
  74. #define udivl __udivsi3 
  75. #define divl __divsi3
  76. #define ddiv __divdf3
  77. #define qmult __qmult
  78. #define dmult __muldf3
  79. #define dneg __negdf2
  80. #define dadd __adddf3
  81. #define dcmp __cmpdf2
  82. #define dtoul __fixunsdfsi
  83. #define dtol __fixdfsi
  84. #define ltod __floatsidf
  85. #endif
  86.  
  87. #ifdef L_umulsi3
  88.  
  89. /*_umulsi3 (a, b) unsigned a, b; { return a * b; } */
  90.  
  91. unsigned long umultl(a,b)
  92. unsigned long a, b;
  93. {
  94.     register unsigned long d7, d6, d5, d4;
  95.     
  96.     d7 = a;
  97.     d6 = b;
  98.     d5 = d6;
  99.     d4 = d6;
  100.     /* without the next line, gcc may core dump. Gcc sometimes
  101.        gets confused if you have too many registers */
  102.     
  103.     &a; &b;
  104.     
  105.     /*printf("a %u b %u\n", a, b);*/
  106.     
  107.     /* low word */
  108.     MUL(d7, d6);
  109.     SWAP(d5);
  110.     MUL(d7, d5);
  111.     SWAP(d7);
  112.     MUL(d7, d4);
  113.     d4 += d5;
  114.     SWAP(d4);
  115.     d4 &= 0xFFFF0000;
  116.     d4 += d6;
  117.     return(d4);
  118. }
  119. #endif
  120.  
  121. #ifdef L_mulsi3
  122. /* _mulsi3 (a, b) int a, b; { return a * b; } */
  123.  
  124. long multl(a, b)
  125. long a, b;
  126. {
  127.     int sign = 0;
  128.     long umultl();
  129.     if ((a ^ b) < 0) sign = 1;
  130.     if (a < 0) a = -a;
  131.     if (b < 0) b = -b;
  132.     /*printf("a is %d b is %d\n", a, b);*/
  133.     if (sign) return(- umultl(a,b));
  134.     else return(umultl(a,b));
  135. }
  136.  
  137. #endif
  138.  
  139. #ifdef L_udivsi3
  140. /*_udivsi3 (a, b) unsigned a, b; { return a / b; } */
  141.  
  142. /*
  143.   this routine based on one in the PD forth package for the sun by Mitch Bradley
  144.   */
  145.  
  146. unsigned udivl(u, v)
  147. register unsigned long u, v;
  148. {
  149.     register unsigned short um, ul, vm, vl;
  150.     unsigned long ans;
  151.     unsigned long u1, v1;
  152.     long i;
  153.     long rem;
  154.     
  155.     if (v == 0) {
  156.     /* should cause an exception condition */
  157.     DIV(u, v);
  158.     fprintf(stderr, "division by zero\n");
  159.     }
  160.     if (v > u) return(0);
  161.     
  162.     ul = u; SWAP(u); um = u;
  163.     vl = v; SWAP(v); vm = v;
  164.     if (vm == 0) {
  165.     u = vl; v = um;
  166.     DIV(u, v);
  167.     /* note -- if you delete the next line, gcc goes into
  168.        an infinite loop */
  169.     if (vm) printf("result is %d\n", v);
  170.     vm = v & 0xFFFF; /* dividend is in low word */
  171.     v &= 0xFFFF0000; /* remainder is in high word */
  172.     v += ul;
  173.     DIV(vl, v);
  174.     v &= 0xFFFF; /* dividend is in low word */
  175.     u = vm;
  176.     SWAP(u);
  177.     return(v + u);
  178.     /*ans = ((um / vl) << 16) + ((((um % vl) << 16) + ul) / vl);
  179.       return(ans);*/
  180.     }
  181.     
  182.     if (vl == 0) return(um / vm);
  183.     SWAP(u); SWAP(v);
  184.     if ( (u >> 3) < v) {
  185.     for(i = 0; u >= v; i++, u -= v);
  186.     /*printf("lt 8\n");*/
  187.     return(i);
  188.     }
  189.     u1 = u; v1 = v;
  190.     
  191.     /* scale divisor */
  192.     v1--;
  193.     for(i = 0; ((unsigned) v1) >= 0xFFFF; v1 >>= 1, i++);
  194.     if (++v1 > 0xFFFF) {
  195.     i++; v1 >>=1;
  196.     }
  197.     u1 >>= i;
  198.     /*printf("i is %d, u1 is %x, v1 is %x\n", i, u1, v1);*/
  199.     ans = u1 / v1;
  200.     rem = u - (ans * v);
  201.     if (rem > v) {ans++; rem -= v; }
  202.     if (rem > v) {printf("oops\n");}
  203.     return(ans);
  204. }
  205. #endif
  206.  
  207. #ifdef L_divsi3
  208.  
  209. long divl(a, b)
  210. long a, b;
  211. {
  212.     int sign = 0;
  213.     if ((a ^ b) < 0) sign = 1;
  214.     if (a < 0) a = -a;
  215.     if (b < 0) b = -b;
  216.     if (sign) return(-udivl(a,b));
  217.     else return(udivl(a,b));
  218. }
  219. #endif
  220.  
  221. #ifdef L_umodsi3
  222. #ifdef __GCC_OLD__
  223. _umodsi3 (a, b)
  224. #else
  225. __umodsi3 (a, b)
  226. #endif
  227. unsigned a, b;
  228. {
  229.     /*return a % b;*/
  230.     return (a - ((a/b)*b));
  231. }
  232. #endif
  233.  
  234. #ifdef L_modsi3
  235. #ifdef __GCC_OLD__
  236. _modsi3 (a, b)
  237. #else
  238. __modsi3 (a, b)
  239. #endif
  240. int a, b;
  241. {
  242.     /*return a % b;*/
  243.     return( a - ((a/b) * b));
  244. }
  245. #endif
  246.  
  247. #ifdef L_lshrsi3
  248. #ifdef __GCC_OLD__
  249. _lshrsi3 (a, b)
  250. #else
  251. __lshrsi3 (a, b)
  252. #endif
  253. unsigned a, b;
  254. {
  255.     return a >> b;
  256. }
  257. #endif
  258.  
  259. #ifdef L_lshlsi3
  260. #ifdef __GCC_OLD__
  261. _lshlsi3 (a, b)
  262. #else
  263. __lshlsi3 (a, b)
  264. #endif
  265. unsigned a, b;
  266. {
  267.     return a << b;
  268. }
  269. #endif
  270.  
  271. #ifdef L_ashrsi3
  272. #ifdef __GCC_OLD__
  273. _ashrsi3 (a, b)
  274. #else
  275. __ashrsi3 (a, b)
  276. #endif
  277. int a, b;
  278. {
  279.     return a >> b;
  280. }
  281. #endif
  282.  
  283. #ifdef L_ashlsi3
  284. #ifdef __GCC_OLD__
  285. _ashlsi3 (a, b)
  286. #else
  287. __ashlsi3 (a, b)
  288. #endif
  289. int a, b;
  290. {
  291.     return a << b;
  292. }
  293. #endif
  294.  
  295.  
  296. /* this divide stuff is hopelessly broken, in addition to being much
  297.    more complicated than in needs to be.  The new version's at the
  298.    end of this file */
  299.  
  300. #if 0
  301. #ifdef L_divdf3
  302.  
  303. /*double _divdf3 (a, b) double a, b; { return a / b; } */
  304.  
  305. double drecip1(f1)
  306. double f1;
  307. {
  308.     struct bitdouble *bdp = &f1;
  309.     unsigned m1, m2;
  310.     
  311.     printf("drecip1(%X)", f1);
  312.     if (bdp->exp == 0 ) return(0L);
  313.     if (bdp->mant1 == 0L) {
  314.     bdp->exp = 0x3ff + 0x3ff - bdp->exp;
  315.     bdp->mant2 = 0L;
  316.     return(f1);
  317.     }
  318.     bdp->exp = 0x3ff + 0x3ff - bdp->exp - 1;
  319.     m1 = (0x00100000 + bdp->mant1) >> 5;
  320.     m2 = (0x80000000 / m1);
  321.     /*printf("m1 %x m2 %x\n", m1, m2);*/
  322.     m2 <<= 5;
  323.     m2 &= 0xFFFFF;
  324.     /*printf("exp %x mant %x\n", bdp->exp, m2);*/
  325.     bdp->mant1 = m2;
  326.     bdp->mant2 = 0L;
  327.     printf("drecip1->%X\n", f1);
  328.     return(f1);
  329. }
  330.  
  331. double drecip(f)
  332. double f;
  333. {
  334.     struct bitdouble *bdp;
  335.     double quotient, remainder;
  336.     
  337.     printf("drecip(%X)", f);
  338.     quotient = drecip1(f);
  339.     remainder = /* 1.0 */ ((double)one) - quotient * f;
  340.     bdp = &remainder;
  341.     for(; bdp->exp > 0x3ca; ) {
  342.     printf("drc: exp=%X ", bdp->exp);
  343.     remainder = /* 1.0 */ ((double)one) - (quotient*f);
  344.     printf("rem=%X ", remainder);
  345.     quotient = quotient + (drecip1(f)*remainder);
  346.     }
  347.     printf("drecip->%X\n", quotient);
  348.     return(quotient);
  349. }
  350.  
  351.  
  352. double ddiv(f1, f2)
  353. double f1, f2;
  354. {
  355.     return(f1 * drecip(f2));
  356. }
  357. #endif
  358. #endif            /* commented out divide routines */
  359.  
  360. #ifdef L_muldf3
  361. /*double _muldf3 (a, b) double a, b; { return a * b; } */
  362.  
  363. #ifdef SLOW_KLUDGEY_QMULT
  364. __qmult(m11, m12, m21, m22, p1, p2)
  365. unsigned long m11, m12, m21, m22, *p1, *p2;
  366. {
  367. /*    register unsigned long d2 = m11; */
  368.     register long d2 = m11;
  369.     register unsigned long d3 = m12, d4 = m21, d5 = m22, d6 =0, d7 = 0;
  370.     int i;
  371.     /* guess what happens if you delete the next line.... */
  372.     /*    &i; */
  373.     for (i = 0; i < 11; i++) ASL2(d2, d3);
  374.     for (i = 0; i < 9; i++) ASL2(d4, d5);
  375.     
  376.     for (i = 0; i < 64; i++) {
  377.     if (d2 < 0) { ADD2(d4, d5, d6, d7);}
  378.     ASL2(d2, d3);
  379.     ASR2(d4, d5);
  380.     }
  381.     d2 = d6;
  382.     d3 = d7;
  383.     for (i = 0; i < 9; i++) ASR2(d2, d3);
  384.     *p1 = d2; *p2 = d3;
  385. }
  386.  
  387. #else
  388. /* new qmult */
  389.  
  390. /* a set of totally local kludges, for summing partial results into
  391.    result vector */
  392. /*
  393. #define XADDL(partial, target_ptr) \
  394.     { register unsigned long temp = *target_ptr; \
  395.     asm volatile("addl %2,%0" : "=d" (temp) : "0" (temp), "g" (partial)); \
  396.     *target_ptr-- = temp; temp = *target_ptr; \
  397.     asm volatile("addxl #0,%0" : "=d" (temp) : "0" (temp)); \
  398.     *target_ptr = temp; }
  399. */
  400.  
  401. static unsigned long constant_zero_kludge = 0;
  402. #define XXADDL(partial, target) \
  403.     { register unsigned long * zero = &constant_zero_kludge + 1; \
  404.       asm volatile("addl  %3,%0@;\
  405.             addxl %1@-,%0@-" \
  406.             : "=a" (target), "=a" (zero) \
  407.             : "0" (target), "g" (partial), "1" (zero)); \
  408.     }
  409.  
  410. /*
  411. #define ADDL(partial, target_ptr) \
  412.     { register unsigned long temp = *target_ptr; \
  413.     asm volatile("addl %2,%0" : "=d" (temp) : "0" (temp), "g" (partial)); \
  414.     *target_ptr-- = temp }
  415.  
  416. #define ADDXL(partial, target_ptr) \
  417.     { register unsigned long temp = *target_ptr; \
  418.     asm volatile("addxl %2,%0" : "=d" (temp) : "0" (temp), "g" (partial)); \
  419.     *target_ptr-- = temp }
  420.     
  421. #define ADDW(partial, target_ptr) \
  422.     { register unsigned short temp = *(unsigned short * )target_ptr; \
  423.     asm volatile("addw %2,%0" : "=d" (temp) : "0" (temp), "g" (partial)); \
  424.     *(unsigned short * )target_ptr-- = temp }
  425.  
  426. #define ADDXW(partial, target_ptr) \
  427.     { register unsigned sort temp = *(unsigned short * )target_ptr; \
  428.     asm volatile("addxw %2,%0" : "=d" (temp) : "0" (temp), "g" (partial)); \
  429.     *(unsigned short * )target_ptr-- = temp }
  430. */    
  431.  
  432. static char component_index[] = 
  433.     {
  434.     3, 3,                /* last ones */
  435.  
  436.     2, 3,                /* next to last x last */
  437.     3, 2,                /* ... */
  438.  
  439.     1, 3,
  440.     2, 2,
  441.     3, 1,
  442.  
  443.     0, 3,
  444.     1, 2,
  445.     2, 1,
  446.     3, 0,
  447.     
  448.     0, 2,
  449.     1, 1,
  450.     2, 0,
  451.     
  452.     0, 1,
  453.     1, 0,
  454.     
  455.     0, 0
  456.     };
  457.  
  458. qmult(m1_hi, m1_lo, m2_hi, m2_lo, result_hi, result_lo)
  459. unsigned long m1_hi, m1_lo, m2_hi, m2_lo, * result_hi, * result_lo;
  460. {
  461.   unsigned short * m1 = (unsigned short * )(&m1_hi);
  462.   unsigned short * m2 = (unsigned short * )(&m2_hi);
  463.   unsigned short result[10];        /* two extra for XADDL */
  464.   register unsigned long mult_register;
  465.   register unsigned long res1, res2, res3;
  466.   long * kludge;
  467.   short i, j;
  468.   char * idx_p = (char * )&component_index;
  469. /*
  470. fprintf(stderr, "  qmult: %08X:%08X * %08X:%08X\n", 
  471.     m1_hi, m1_lo, m2_hi, m2_lo);
  472. */
  473.   for (i = 0 ; i < 10 ; i++) result[i] = 0;
  474.  
  475. /* walk thru 'vectors' of mantissa pieces, doing unsigned multiplies
  476.    and summing results into result vector.  Note that the order is 
  477.    chosen to guarantee that no more adds than generated by XADDL are
  478.    needed.  */
  479.  
  480.   for ( ; ; )
  481.     {
  482.     i = *idx_p++;
  483.     j = *idx_p++;
  484.     mult_register = m1[i]; 
  485.     MUL(m2[j], mult_register);
  486.     kludge = (long * )(&result[i + j + 2]);
  487.     XXADDL(mult_register, kludge);
  488. /*
  489. fprintf(stderr, "  qmult: %d[%04X] * %d[%04X] -> %08X\n",
  490.   i, m1[i], j, m2[j], mult_register);
  491. fprintf(stderr, "       %04X %04X %04X %04X %04X %04X %04X %04X\n",
  492.   result[2], result[3], result[4], result[5], 
  493.   result[6], result[7], result[8], result[9]);
  494. */
  495.     if ((i == 0) && (j == 0))
  496.         break;
  497.     }
  498.  
  499.   /* get the result shifted to the right place.  Unfortunately, we need
  500.      the 53 bits that's 22 bits down in the result vec.  sigh */
  501.   kludge = (long * )(&result[2]);
  502.   res1 = *kludge++;
  503.   res2 = *kludge++;
  504.   res3 = *kludge;
  505.   for (i = 0 ; i < 12 ; i++)
  506.     ASL3(res1, res2, res3);
  507.   /* now put the result back */
  508.   *result_hi = res1;
  509.   *result_lo = res2;
  510. }
  511. #endif
  512.  
  513.  
  514. double dmult(f1, f2)
  515. double f1, f2;
  516. {
  517.     register long d2;
  518.     register unsigned d3, d4, d5, d6, d7;
  519.     unsigned long p1, p2;
  520.     
  521.     struct bitdouble
  522.     *bdp1 = (struct bitdouble *)&f1,
  523.     *bdp2 = (struct bitdouble *)&f2;
  524.     int exp1, exp2, i;
  525.     
  526.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  527.     /* check for zero */
  528.     if (! exp1) return(0.0);
  529.     if (! exp2) return(0.0);
  530.     d2 = 0x00100000 + bdp1->mant1;
  531.     d3 = bdp1->mant2;
  532.     d4 = 0x00100000 + bdp2->mant1;
  533.     d5 = bdp2->mant2;
  534.     __qmult(d2, d3, d4, d5, &p1, &p2);
  535.     d6 = p1; d7 = p2;
  536.     
  537.     if (d6 & 0x00200000) {
  538.     ASR2(d6, d7);
  539.     exp1++;
  540.     }
  541.     
  542.     if (bdp1->sign == bdp2->sign) bdp1->sign = 0;
  543.     else bdp1->sign = 1;
  544.     bdp1->exp = exp1 + exp2 - 0x3ff;
  545.     bdp1->mant1 = d6;
  546.     bdp1->mant2 = d7;
  547.     return(f1);
  548. }
  549. #endif
  550.  
  551. #ifdef L_negdf2
  552. /*double _negdf2 (a) double a; { return -a; } */
  553.  
  554. double dneg(num)
  555. double num;
  556. {
  557.     long *i = (long *)#
  558.     if (*i)            /* only negate if non-zero */
  559.         *i ^= 0x80000000;
  560.     return(num);
  561. }
  562. #endif
  563.  
  564. #ifdef L_adddf3
  565. /*double _adddf3 (a, b) double a, b; { return a + b; } */
  566.  
  567. double dadd(f1, f2)
  568. double f1, f2;
  569. {
  570.     
  571.     register long d4, d5, d6, d7;
  572.     struct bitdouble
  573.     *bdp1 = (struct bitdouble *)&f1,
  574.         *bdp2 = (struct bitdouble *)&f2;
  575.     short exp1, exp2, sign1, sign2, howmuch, temp;
  576.     
  577.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  578.     
  579.     /* check for zero */
  580.     if (! exp1) return(f2); if (! exp2) return(f1);
  581.     
  582.     /* align them */
  583.     if (exp1 < exp2) {
  584.     bdp1 = (struct bitdouble *)&f2; bdp2 = (struct bitdouble *)&f1;
  585.     exp1 = bdp1->exp; exp2 = bdp2->exp;
  586.     }
  587.     howmuch = exp1 - exp2;
  588.     if (howmuch > 53) return(f1);
  589.     
  590.     d7 = bdp2->mant1 + 0x00100000;
  591.     d6 = bdp2->mant2;
  592.     
  593.     d5 = bdp1->mant1 + 0x00100000;
  594.     d4 = bdp1->mant2;
  595.     
  596.     for (temp = 0; temp < howmuch; temp++) ASR2(d7, d6);
  597.     
  598.     /* take care of negative signs */
  599.     if (bdp1->sign) 
  600.     {
  601.     NEG(d5, d4);
  602.     }
  603.     if (bdp2->sign)
  604.     {
  605.     NEG(d7, d6);
  606.     }
  607.     
  608.     ADD2(d7, d6, d5, d4);
  609.     bdp1 = (struct bitdouble *)&f1;
  610.     
  611.     if (d5 < 0) {
  612.     NEG(d5, d4);
  613.     bdp1->sign = 1;
  614.     } else bdp1->sign = 0;
  615.     
  616.     if (d5 == 0 && d4 == 0) return(0.0);
  617.     
  618.     for(howmuch = 0; d5 >= 0; howmuch++) ASL2(d5, d4);
  619.     
  620.     ASL2(d5, d4);
  621.     for (temp = 0; temp < 12; temp++) ASR2(d5, d4);
  622.     bdp1->mant1 = d5;
  623.     bdp1->mant2 = d4;
  624.     bdp1->exp = exp1 + 11 - howmuch;
  625.     return(f1); 
  626. }
  627.  
  628. #endif
  629.  
  630. #ifdef L_subdf3
  631. double
  632. #ifdef __GCC_OLD__
  633.     _subdf3 (a, b)
  634. #else
  635.     __subdf3 (a, b)
  636. #endif
  637. double a, b;
  638. {
  639.     return a+(-b);
  640. }
  641. #endif
  642.  
  643. #ifdef L_cmpdf2
  644. /*
  645.   int _cmpdf2 (a, b) double a, b; { if (a > b) return 1;
  646.   else if (a < b) return -1; return 0; } 
  647.   */
  648.  
  649. int dcmp(f1, f2)
  650. double f1, f2;
  651. {
  652.     struct bitdouble *bdp1, *bdp2;
  653.     unsigned int s1, s2;
  654.     bdp1 = (struct bitdouble *)&f1;
  655.     bdp2 = (struct bitdouble *)&f2;
  656.     s1 = bdp1->sign;
  657.     s2 = bdp2->sign;
  658.     if (s1 > s2) return(-1);
  659.     if (s1 < s2) return(1);
  660.     /* if sign of both negative, switch them */
  661.     if (s1 != 0) {
  662.     bdp1 = (struct bitdouble *)&f2;
  663.     bdp2 = (struct bitdouble *)&f1;
  664.     }
  665.     s1 = bdp1->exp;
  666.     s2 = bdp2->exp;
  667.     if (s1 > s2) return(1);
  668.     if (s1 < s2) return(-1);
  669.     /* same exponent -- have to compare mantissa */
  670.     s1 = bdp1->mant1;
  671.     s2 = bdp2->mant1;
  672.     if (s1 > s2) return(1);
  673.     if (s1 < s2) return(-1);
  674.     s1 = bdp1->mant2;
  675.     s2 = bdp2->mant2;
  676.     if (s1 > s2) return(1);
  677.     if (s1 < s2) return(-1);
  678.     return(0); /* the same! */
  679. }
  680. #endif
  681.  
  682. #ifdef L_fixunsdfsi
  683. /*_fixunsdfsi (a) double a; { return (unsigned int) a; } */
  684.  
  685. /* #ifdef L_fixdfsi _fixdfsi (a) double a; { return (int) a; } #endif */
  686.  
  687. unsigned long dtoul(f)
  688. double f;
  689. {
  690.     struct bitdouble *bdp;
  691.     int si, ex, mant1, mant2;
  692.     bdp = (struct bitdouble *)&f;
  693.     si = bdp->sign;
  694.     ex = bdp->exp;
  695.     mant1 = bdp->mant1 + 0x00100000;
  696.     mant2 = bdp->mant2;
  697.     
  698.     /* zero value */
  699.     if (ex == 0) return(0L);
  700.     /* less than 1 */
  701.     if (ex < 0x3ff) return(0L);
  702.     /* less than 0 */
  703.     if (si ) return(0L);
  704.     mant1 <<= 10;
  705.     mant1 += (mant2 >> 22);
  706.     mant1 >>= 30 + (0x3ff - ex);
  707.     return(mant1);
  708. }
  709.  
  710. long dtol(f)
  711. double f;
  712. {
  713.     struct bitdouble *bdp = (struct bitdouble *)&f;
  714.     
  715.     if (bdp->sign) {
  716.     bdp->sign = 0;
  717.     return( - dtoul(f));
  718.     }
  719.     return(dtoul(f));
  720. }
  721. #endif
  722.  
  723. #ifdef L_fixunsdfdi
  724.  
  725. /*double _fixunsdfdi (a) double a; { union double_di u;
  726.   u.i[LOW] = (unsigned int) a; u.i[HIGH] = 0; return u.d; } */
  727. #endif
  728.  
  729.  
  730. #ifdef L_fixdfdi
  731. double
  732. #ifdef __GCC_OLD__
  733.     _fixdfdi (a)
  734. #else
  735.     __fixdfdi (a)
  736. #endif
  737. double a;
  738. {
  739.     union double_di u;
  740.     u.i[LOW] = (int) a;
  741.     u.i[HIGH] = (int) a < 0 ? -1 : 0;
  742.     return u.d;
  743. }
  744. #endif
  745.  
  746. #ifdef L_floatsidf
  747. /* double _floatsidf (a) int a; { return (double) a; } */
  748.  
  749. double ltod(i)
  750. long i;
  751. {
  752.     int expo, shift;
  753.     double retval;
  754.     struct bitdouble *bdp = (struct bitdouble *)&retval;
  755.     if (i == 0) {
  756.     long *t = (long *)&retval;
  757.     t[0] = 0L;
  758.     t[1] = 0L;
  759.     return(retval);
  760.     }
  761.     if (i < 0) {
  762.     bdp->sign = 1;
  763.     i = -i;
  764.     } else bdp->sign = 0;
  765.     shift = i;
  766.     for (expo = 0x3ff + 31 ; shift > 0; expo--, shift <<= 1);
  767.     shift <<= 1;
  768.     bdp->exp = expo;
  769.     bdp->mant1 = shift >> 12;
  770.     bdp->mant2 = shift << 20;
  771.     return(retval);
  772. }
  773.  
  774. #endif
  775.  
  776. #ifdef L_floatdidf
  777. /* ok as is -- will call other routines */
  778. double
  779. #ifdef __GCC_OLD__
  780.     _floatdidf (u)
  781. #else
  782.     __floatdidf (u)
  783. #endif
  784. union double_di u;
  785. {
  786.     register double hi
  787.     = ((double) u.i[HIGH]) * (double) 0x10000 * (double) 0x10000;
  788.     register double low = (unsigned int) u.i[LOW];
  789.     return hi + low;
  790. }
  791. #endif
  792.  
  793. #ifdef L_addsf3
  794. SFVALUE
  795.     _addsf3 (a, b)
  796. union flt_or_int a, b;
  797. {
  798.     union flt_or_int intify; return INTIFY ((double) a.f + (double) b.f);
  799. }
  800. #endif
  801.  
  802. #ifdef L_negsf2
  803. SFVALUE
  804.     _negsf2 (a)
  805. union flt_or_int a;
  806. {
  807.     union flt_or_int intify;
  808.     return INTIFY (-((double) a.f));
  809. }
  810. #endif
  811.  
  812. #ifdef L_subsf3
  813. SFVALUE
  814.     _subsf3 (a, b)
  815. union flt_or_int a, b;
  816. {
  817.     union flt_or_int intify;
  818.     return INTIFY (((double) a.f - (double) b.f));
  819. }
  820. #endif
  821.  
  822. #ifdef L_cmpsf2
  823. SFVALUE
  824.     _cmpsf2 (a, b)
  825. union flt_or_int a, b;
  826. {
  827.     union flt_or_int intify;
  828.     double a1, b1;
  829.     a1 = a.f; b1 = b.f;
  830.     if ( a1 > b1)
  831.     return 1;
  832.     else if (a1 < b1)
  833.     return -1;
  834.     return 0;
  835. }
  836. #endif
  837.  
  838. #ifdef L_mulsf3
  839. SFVALUE
  840.     _mulsf3 (a, b)
  841. union flt_or_int a, b;
  842. {
  843.     union flt_or_int intify;
  844.     return INTIFY (((double) a.f * (double) b.f));
  845. }
  846. #endif
  847.  
  848. #ifdef L_divsf3
  849. SFVALUE
  850.     _divsf3 (a, b)
  851. union flt_or_int a, b;
  852. {
  853.     union flt_or_int intify;
  854.     return INTIFY (((double) a.f / (double) b.f));
  855. }
  856. #endif
  857.  
  858. #ifdef L_truncdfsf2
  859. float __dtof(d)
  860. double d;
  861. {
  862.     struct bitdouble *bdp = (struct bitdouble *)&d;
  863.     float retval;
  864.     struct bitfloat *bfp = (struct bitfloat *)&retval;
  865.     int tempval;
  866.     
  867.     bfp->sign = bdp->sign;
  868.     if (bdp->exp == 0) return ((float) 0.0);
  869.     bfp->exp = bdp->exp - 0x400 + 0x80;
  870.     tempval = (bdp->mant1 << 4 ) + ((0xF0000000 & bdp->mant2) >> 28);
  871.     /* round */
  872.     tempval++;
  873.     if (tempval == 0x01000000) bfp->exp++;
  874.     bfp->mant = tempval >> 1;
  875.     return(retval);
  876. }
  877.  
  878. SFVALUE
  879. #ifdef __GCC_OLD__
  880.     _truncdfsf2 (a)
  881. #else
  882.     __truncdfsf2 (a)
  883. #endif
  884. double a;
  885. {
  886.     union flt_or_int intify;
  887.     return INTIFY (__dtof(a));
  888. }
  889. #endif
  890.  
  891. #ifdef L_extendsfdf2
  892. double __ftod(f)
  893. union flt_or_int f;
  894. {
  895.     double retval;
  896.     struct bitfloat *bfp = (struct bitfloat *)&f.f;
  897.     struct bitdouble *bdp = (struct bitdouble *)&retval;
  898.     if (bfp->exp == 0) return(0.0);
  899.     bdp->sign = bfp->sign;
  900.     bdp->exp = 0x400 - 0x80 + bfp->exp;
  901.     bdp->mant1 = bfp->mant >> 3;
  902.     bdp->mant2 = (bfp->mant & 0x7) << 29;
  903.     /*printf("returning %f from extendsfdf2\n", retval);*/
  904.     return(retval);
  905. }
  906.  
  907. double
  908. #ifdef __GCC_OLD__
  909.     _extendsfdf2 (a)
  910. #else
  911.     __extendsfdf2 (a)
  912. #endif
  913. union flt_or_int a;
  914. {
  915.     union flt_or_int intify;
  916.     double retval;
  917.     return (__ftod(a));
  918. }
  919. #endif
  920.  
  921. #ifdef L_divdf3
  922.  
  923. /* new double-float divide routine, by jrd */
  924. /* thanks jrd !! */
  925.  
  926. #ifdef __GCC_OLD__
  927. double _divdf3(num, denom)
  928. #else
  929. double __divdf3(num, denom)
  930. #endif
  931. double num, denom;
  932. {
  933.     double local_num = num;
  934.     double local_denom = denom;
  935.     struct bitdouble * num_bdp = (struct bitdouble *)(&local_num);
  936.     struct bitdouble * den_bdp = (struct bitdouble *)(&local_denom);
  937.     short num_exp = num_bdp->exp,
  938.           den_exp = den_bdp->exp;
  939.     short result_sign = 0;
  940.     /*  register */ unsigned long num_m1, num_m2, num_m3, num_m4;
  941.     register unsigned long den_m1, den_m2, den_m3, den_m4;
  942.     unsigned long result_mant[3];
  943.     unsigned long result_mask;
  944.     short result_idx;
  945.     
  946.     if ((num_exp == 0) || (den_exp == 0))        /* zzz should really cause trap */
  947.       return(0.0);
  948.     
  949.     /* deal with signs */
  950.     result_sign = result_sign + num_bdp->sign - den_bdp->sign;
  951.     
  952.     /* unpack the numbers */
  953.     num_m1 = num_bdp->mant1 | 0x00100000;        /* hidden bit */
  954.     num_m2 = num_bdp->mant2;
  955.     num_m4 = num_m3 = 0;
  956.     den_m1 = den_bdp->mant1 | 0x00100000;        /* hidden bit */
  957.     den_m2 = den_bdp->mant2;
  958.     den_m4 = den_m3 = 0;
  959.     
  960. #if 0
  961.     /* buy a little extra accuracy by shifting num and denum up 10 bits */
  962.     for (result_idx /* loop counter */ = 0 ; result_idx < 10 ; result_idx++)
  963.     {
  964.     ASL3(num_m1, num_m2, num_m3);
  965.     ASL3(den_m1, den_m2, den_m3);
  966.     }
  967. #endif
  968.     
  969.     /* hot wire result mask and index */
  970.     result_mask = 0x00100000;            /* start at top mantissa bit */
  971.     result_idx = 0;                /* of first word */
  972.     result_mant[0] = 0;
  973.     result_mant[1] = 0;
  974.     result_mant[2] = 0;
  975.     
  976.     /* if denom is greater than num here, shift denom right one and dec num expt */
  977.     if (den_m1 < num_m1)
  978.       goto kludge1;                /* C is assembly language,
  979.                            remember? */
  980.     if (den_m1 > num_m1)
  981.       goto kludge0;
  982.     if (den_m2 <= num_m2)                /* first word eq, try 2nd */
  983.       goto kludge1;
  984.     
  985.   kludge0:
  986.     
  987.     num_exp--;
  988.     ASR4(den_m1, den_m2, den_m3, den_m4);
  989.     
  990.   kludge1:
  991.     
  992.     for ( ; !((result_idx == 2) && (result_mask & 0x40000000)) ; )
  993.     {
  994.     /* if num >= den, subtract den from num and set bit in result */
  995.     if (num_m1 > den_m1) goto kludge2;
  996.     if (num_m1 < den_m1) goto kludge3;
  997.     if (num_m2 > den_m2) goto kludge2;
  998.     if (num_m2 < den_m2) goto kludge3;
  999.     if (num_m3 > den_m3) goto kludge2;
  1000.     if (num_m3 < den_m3) goto kludge3;
  1001.     if (num_m4 < den_m4) goto kludge3;
  1002.     
  1003.       kludge2:
  1004.     result_mant[result_idx] |= result_mask;
  1005.     SUB4(den_m1, den_m2, den_m3, den_m4, num_m1, num_m2, num_m3, num_m4);
  1006.       kludge3:
  1007.     ASR4(den_m1, den_m2, den_m3, den_m4);
  1008.     result_mask >>= 1;
  1009.     if (result_mask == 0)        /* next word? */
  1010.     {
  1011.         result_mask = 0x80000000;
  1012.         result_idx++;
  1013.     }
  1014.     }
  1015.  
  1016. /* stick in last half-bit if necessary */
  1017.   if (result_mant[2])
  1018.     {
  1019.     den_m1 = 0;        /* handy register */
  1020.     den_m2 = 1;
  1021. /* busted
  1022.     ADD2(den_m1, den_m2, result_mant[0], result_mant[1]);
  1023. */
  1024.     result_mant[1]++;
  1025.     if (result_mant[1] == 0)
  1026.         result_mant[0]++;
  1027.  
  1028.     if (result_mant[0] & 0x00200000)    /* overflow? */
  1029.         {
  1030.         num_exp--;
  1031.         }
  1032.     }
  1033.   /* compute the resultant exponent */
  1034.   num_exp = num_exp - den_exp + 0x3FF;
  1035.  
  1036.   /* reconstruct the result in local_num */
  1037.   num_bdp->sign = result_sign;
  1038.   num_bdp->exp = num_exp;
  1039.   num_bdp->mant1 = result_mant[0] & 0xFFFFF;
  1040.   num_bdp->mant2 = result_mant[1];
  1041.   
  1042.   /* done! */
  1043.   return(local_num);
  1044. }
  1045. #endif
  1046.