home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / C / SCBENCH.ZIP / SCFM87.C < prev    next >
Encoding:
Text File  |  1980-01-01  |  10.7 KB  |  609 lines
  1. /*
  2. ** BYTE Small C Floating Point Coprocessor Library
  3. ** BYTE Small C based on Small C by J.E. Hendrix
  4. ** 8087 version for 8088/8086
  5. ** Version 1, Feb. 17, 1988
  6. **
  7. ** Conventions:
  8. ** All variables are manipulated OUTSIDE the coprocessor in
  9. ** 8-byte form.
  10. */
  11.  
  12. int    fsword;        /* Storage for status word */
  13. int    fsword1;    /* Extra storage for same */
  14. /*
  15. ** Initialize coprocessor
  16. ** Call this at the beginning and end of the program
  17. */
  18. finit() {
  19. #asm
  20.     FINIT
  21. #endasm
  22. }
  23.  
  24. /*
  25. ** fint2float(ptr)
  26. ** Converts integer pointed to by (ptr) into 8-byte floating point
  27. ** form and stores it back into (ptr);
  28. */
  29. fint2float(ptr) int ptr[]; {
  30. #asm
  31.     MOV    BP,SP
  32.     MOV    SI,2[BP]    ;Get pointer to 8-byte block
  33.     FILD    QWORD PTR [SI]    ;Load the block in
  34.     FSTP    QWORD PTR [SI]    ;Put it back out
  35. #endasm
  36. }
  37.  
  38. /*
  39. ** ffloat2int(ptr)
  40. ** Converts 8-byte floating point number at ptr back to an integer.
  41. */
  42. ffloat2int(ptr) int ptr[];
  43. {
  44. #asm
  45.     MOV    BP,SP
  46.     MOV    SI,2[BP]    ;Get pointer to 8-byte block
  47.     FLD    QWORD PTR [SI]    ;Load it in
  48.     FISTP    QWORD PTR [SI]    ;Save it back
  49. #endasm
  50. }
  51.  
  52. /*
  53. ** f2add(fac1,fac2,dest)
  54. ** Add fac1 to fac2, store the result in dest
  55. */
  56. f2add(fac1,fac2,dest) int fac1[],fac2[],dest[];
  57. {
  58. #asm
  59.     MOV    BP,SP
  60.     MOV    SI,6[BP]    ;Get first factor
  61.     FLD    QWORD PTR [SI]    ;Push it
  62.     MOV    SI,4[BP]    ;Get second factor
  63.     FADD    QWORD PTR [SI]    ;Add it
  64.     MOV    SI,2[BP]    ;Get destination
  65.     FSTP    QWORD PTR [SI]    ;Store
  66. #endasm
  67. }
  68.  
  69. /*
  70. ** f2sub(fac1,fac2,dest)
  71. ** Subtract fac2 from fac1, store result in dest
  72. */
  73. f2sub(fac1,fac2,dest) int fac1[],fac2[],dest[];
  74. {
  75. #asm
  76.     MOV    BP,SP
  77.     MOV    SI,6[BP]    ;Get fac1
  78.     FLD    QWORD PTR [SI]    ;Push it
  79.     MOV    SI,4[BP]    ;Get fac2
  80.     FSUB    QWORD PTR [SI]    ;Subtract
  81.     MOV    SI,2[BP]    ;Destination
  82.     FSTP    QWORD PTR [SI]    ;Store
  83. #endasm
  84. }
  85.  
  86. /*
  87. ** f2mult(fac1,fac2,dest)
  88. ** Multiply fac1 by fac2, result in dest
  89. */
  90. f2mult(fac1,fac2,dest) int fac1[],fac2[],dest[];
  91. {
  92. #asm
  93.     MOV    BP,SP
  94.     MOV    SI,6[BP]    ;get fac1
  95.     FLD    QWORD PTR [SI]    ;Push it
  96.     MOV    SI,4[BP]    ;Get fac2
  97.     FLD    QWORD PTR [SI]
  98.     FMUL
  99.     MOV    SI,2[BP]    ;Destination pointer
  100.     FSTP    QWORD PTR [SI]    ;Store
  101. #endasm
  102. }
  103.  
  104. /*
  105. ** f2div(fac1,fac2,dest)
  106. ** Divide fac1 by fac2, result in dest
  107. */
  108. f2div(fac1,fac2,dest) int fac1[],fac2[],dest[];
  109. {
  110. #asm
  111.     MOV    BP,SP
  112.     MOV    SI,6[BP]    ;Get fac1
  113.     FLD    QWORD PTR [SI]    ;Push it
  114.     MOV    SI,4[BP]    ;Get fac2
  115.     FDIV    QWORD PTR [SI]    ;Divide
  116.     MOV    SI,2[BP]    ;Destination
  117.     FSTP    QWORD PTR [SI]    ;Store it
  118. #endasm
  119. }
  120.  
  121. /*
  122. ** fload(ptr)
  123. ** Load floating point number given by ptr
  124. */
  125. fload(ptr) int ptr[];
  126. {
  127. #asm
  128.     MOV    BP,SP
  129.     MOV    SI,2[BP]
  130.     FLD    QWORD PTR [SI]
  131. #endasm
  132. }
  133.  
  134. /*
  135. ** fstore(ptr)
  136. ** Store floating point number given by ptr - pops floating point stack
  137. */
  138. fstore(ptr) int ptr[];
  139. {
  140. #asm
  141.     MOV    BP,SP
  142.     MOV    SI,2[BP]
  143.     FSTP    QWORD PTR [SI]
  144. #endasm
  145. }
  146.  
  147. /*
  148. ** fadd(ptr)
  149. ** Add contents of ptr to top of floating point stack -- result on stack
  150. */
  151. fadd(ptr) int ptr[];
  152. {
  153. #asm
  154.     MOV    BP,SP
  155.     MOV    SI,2[BP]
  156.     FADD    QWORD PTR [SI]
  157. #endasm
  158. }
  159.  
  160. /*
  161. ** fsub(ptr)
  162. ** Subtract contents of ptr from top of stack -- result on top of stack
  163. */
  164. fsub(ptr) int ptr[];
  165. {
  166. #asm
  167.     MOV    BP,SP
  168.     MOV    SI,2[BP]
  169.     FSUB    QWORD PTR [SI]
  170. #endasm
  171. }
  172.  
  173. /*
  174. ** fmult(ptr)
  175. ** Multiply contents of ptr by top of stack -- result on top of stack
  176. */
  177. fmult(ptr) int ptr[];
  178. {
  179. #asm
  180.     MOV    BP,SP
  181.     MOV    SI,2[BP]
  182.     FMUL    QWORD PTR [SI]
  183. #endasm
  184. }
  185.  
  186. /*
  187. ** fdiv(ptr)
  188. ** Divide top of stack by contents of ptr -- result on top of stack
  189. */
  190. fdiv(ptr) int ptr[];
  191. {
  192. #asm
  193.     MOV    BP,SP
  194.     MOV    SI,2[BP]
  195.     FDIV    QWORD PTR [SI]
  196. #endasm
  197. }
  198.  
  199. /*
  200. ** fabs()
  201. ** Set top of stack to its absolute value
  202. */
  203. fabs()
  204. {
  205. #asm
  206.     FABS
  207. #endasm
  208. }
  209.  
  210. /*
  211. ** fsqrt()
  212. ** Take square root of top of floating point stack
  213. */
  214. fsqrt()
  215. {
  216. #asm
  217.         FSQRT
  218. #endasm
  219. }
  220.  
  221. /*
  222. ** fcompz(ptr)
  223. ** Compares number at ptr with zero
  224. ** Returns: -1 if less, 0 if equal, 1 if greater
  225. */
  226. fcompz(ptr) int ptr[];
  227. {
  228. #asm
  229.     FLDZ            ;get zero on stack
  230.     MOV    BP,SP
  231.     MOV    SI,2[BP]
  232.     FLD    QWORD PTR [SI]    ;Get the number
  233.     FCOMPP            ;Compare
  234.     FSTSW    _fsword        ;Store the control word
  235.     FWAIT            ;Make sure the coprocessor is done
  236.     MOV    AX,_fsword    ;Get the control word
  237.     XOR    BX,BX
  238.     SAHF            ;Move into flags
  239.     JZ    FCZX
  240.     JNC    FCZ1
  241.     DEC    BX
  242.     JMP    FCZX
  243. FCZ1:    INC    BX
  244. FCZX:
  245.  
  246. #endasm
  247. }
  248. /*
  249. ** fcomp(ptr1,ptr2)
  250. ** ptr1 and ptr2 point to floating-point numbers.
  251. ** This function returns:
  252. ** 1 if ptr1>ptr2; 0 if ptr1=ptr2; -1 if ptr1<ptr2
  253. */
  254. fcomp(ptr1,ptr2) int ptr1[],ptr2[];
  255. {
  256. #asm
  257.     MOV    BP,SP
  258.     MOV    SI,2[BP]    ;Get ptr2
  259.     FLD    QWORD PTR [SI]    ;Load it in
  260.  
  261.     MOV    SI,4[BP]
  262.     FLD    QWORD PTR [SI]    ;Get ptr1
  263.     FCOMPP            ;Compare
  264.     FSTSW    _fsword        ;Store the control word
  265.     FWAIT            ;Make sure the coprocessor is done
  266.     MOV    AX,_fsword    ;Get the control word
  267.     XOR    BX,BX
  268.     SAHF            ;Move into flags
  269.     JZ    FCMX
  270.     JNC    FCM1
  271.     DEC    BX
  272.     JMP    FCMX
  273. FCM1:    INC    BX
  274. FCMX:
  275. #endasm
  276. }
  277.  
  278. /*
  279. ** ffact(n) int n;
  280. ** Calculate n!.  n! is left on top of stack.
  281. */
  282. ffact(n) int n;
  283. {
  284. #asm
  285.     MOV    BP,SP        ;get value
  286.     MOV    AX,2[BP]
  287.     FLD1
  288.     FLD1
  289. FFCT0:    DEC    AX        ;done?
  290.     JZ    FFCT1
  291.     FLD1
  292.     FADD    ST,ST(2)    ;bump
  293.     FST    ST(2)
  294.     FMUL            ;Multiply
  295.     JMP    SHORT FFCT0
  296. FFCT1:    FFREE    ST(1)        ;Clear n
  297. #endasm
  298. }
  299.  
  300. /*
  301. ** fsin(ptr,ptr1)
  302. ** Calculate sine of quantity at ptr...returns sine in ptr1
  303. */
  304. fsin(ptr,ptr1) int ptr[],ptr1[];
  305. {
  306. #asm
  307.     MOV    BP,SP
  308.     MOV    SI,4[BP]    ;Get pointer
  309.     FLD    QWORD PTR [SI]    ;Get argument
  310. ;See if arg is 0 -- if so, we can skip a lot
  311.     FTST
  312.     FSTSW    _fsword
  313.     FWAIT
  314.     MOV    AX,_fsword
  315.     SAHF
  316.     JZ    FSIN1
  317.     FFREE    ST(0)
  318.     FLD1            ;Get a 1
  319.     FLD    ST(0)        ;Dup it
  320.     FADD            ;2 now at stack top
  321.     FLD    ST(0)        ;Dup it
  322.     FLD    QWORD PTR [SI]    ;Get argument
  323.     FDIVR            ;arg/2
  324.     FPTAN            ;Get partial tangent
  325.     FDIVR            ;Calc y/x
  326.     FST    ST(2)
  327.     FMUL            ;2*(y/x) on top
  328.     FXCH
  329.     FMUL    ST(0),ST    ;(y/x)^2
  330.     FLD1            ;Get 1 again
  331.     FADD            ;(1+(y/x)^2)
  332.     FDIV
  333. FSIN1:    MOV    SI,2[BP]
  334.     FSTP    QWORD PTR [SI]
  335. #endasm
  336. }
  337.  
  338. /*
  339. ** fxtoy(x,y,z)
  340. ** Given x and y are pointers to floating point numbers,
  341. ** calculates x^y and returns the result in z.
  342. ** Not that this routine is only good for returning roots
  343. ** of things.
  344. */
  345. fxtoy(x,y,z) int x[],y[],z[];
  346. {
  347. #asm
  348.     MOV    BP,SP
  349.     MOV    SI,4[BP]    ;Get y
  350.     FLD    QWORD PTR [SI]    ;y in st(0)
  351.     MOV    SI,6[BP]
  352.     FLD    QWORD PTR [SI]    ;x on top
  353.     FYL2X            ;y*log2(x)
  354.     CALL    _ftwo2x        ;2^(y*log2(x))
  355.     MOV    SI,2[BP]    ;Get z address
  356.     FSTP    QWORD PTR [SI]    ;Store
  357. #endasm
  358. }
  359.  
  360. /*
  361. ** fex(ptr,ptr1)
  362. ** Raises e^x, where x is in ptr.  Places result in ptr1.
  363. **/
  364. fex(ptr,ptr1) int ptr[],ptr1[];
  365. {
  366. #asm
  367.     MOV    BP,SP
  368.     MOV    SI,4[BP]
  369.     FLD    QWORD PTR [SI]    ;x in st(0)
  370.     FLDL2E            ;log2(e)
  371.     FMUL            ;x*log2(e)
  372.     CALL    _ftwo2x        ;2^(x*log2(e))
  373.     MOV    SI,2[BP]    ;Get destination
  374.     FSTP    QWORD PTR [SI]    ;Store
  375. #endasm
  376. }
  377.  
  378. /*
  379. ** fround(ptr)
  380. ** Assume ptr points to quadword floating point.
  381. ** Rounds contents of ptr to integer and stores result in ptr
  382. */
  383. fround(ptr) int ptr[];
  384. {
  385. #asm
  386.     MOV    BP,SP
  387.     MOV    SI,2[BP]    ;Pointer to si
  388.     FLD    QWORD PTR [SI]
  389.     FRNDINT
  390.     FSTP    QWORD PTR [SI]
  391. #endasm
  392. }
  393.  
  394. /*
  395. ** fstoreb(ptr)
  396. ** Store top of stack as BCD into ptr
  397. */
  398. fstoreb(ptr) int ptr[];
  399. {
  400. #asm
  401.     MOV    BP,SP
  402.     MOV    SI,2[BP]    ;Get pointer
  403.     FBSTP    [SI]        ;Store it
  404.     FWAIT
  405. #endasm
  406. }
  407.  
  408. /*
  409. ** ftwo2x()
  410. ** Raises two to the power of x.
  411. ** x is on top of the stack.
  412. ** Leaves result on top of stack.
  413. */
  414. ftwo2x()
  415. {
  416. #asm
  417. ;First save existing control word
  418.     FNSTCW    _fsword
  419.     MOV    AX,_fsword
  420.     AND    AX,0F3FFH    ;Round toward nearest even
  421.     MOV    _fsword1,AX
  422.     FLDCW    _fsword1    ;New control word
  423. ;Now break exponent into integer and fractional part
  424.     FLD    ST(0)        ;Dup x
  425.     FRNDINT            ;Int part in ST(0)
  426.     FSUB    ST(1),ST    ;Fractional part in ST(1)
  427.     FXCH            ;Switch - int in st(1),
  428.                 ;frac. in st(0)
  429.     FTST            ;frac. < 0 ?
  430.     FSTSW    _fsword1
  431.     FWAIT
  432.     F2XM1            ;2^(f)-1
  433. ;See if fractional part is negative
  434.     MOV    AX,_fsword1
  435.     SAHF
  436.     JNC    T2X1
  437. ;Fractional part is neg.
  438.     FCHS            ;Make it positive
  439.     FLD    ST(0)        ;dup
  440.     FLD1            ;1 on top
  441.     FADD            ;2^(-f) on top
  442.     FDIV            ;(2^(-f)-1)/2^(-f)
  443.     FCHS            ;2^f-1
  444. T2X1:    FLD1
  445.     FADD            ;2^f
  446.     FSCALE            ;2^x
  447.     FFREE    ST(1)        ;Clear out i
  448. ;Restore control word
  449.     FLDCW    _fsword
  450. #endasm
  451. }
  452.  
  453. /*
  454. ** fconst(val) int val;
  455. ** Load top of FPU stack with constant based on val:
  456. ** val     constant
  457. ** 0        0.0
  458. ** 1        1.0
  459. ** 2        pi
  460. ** 3        log2(e)
  461. ** 4        log2(10)
  462. ** 5        log10(2)
  463. ** 6        loge(2)
  464. */
  465. fconst(val) int val;
  466. {
  467.     switch(val) {
  468.      default:
  469.      case 0:
  470. #asm
  471.       FLDZ
  472. #endasm
  473.        break;
  474.  
  475.      case 1:
  476. #asm
  477.       FLD1
  478. #endasm
  479.        break;
  480.  
  481.      case 2:
  482. #asm
  483.       FLDPI
  484. #endasm
  485.        break;
  486.  
  487.      case 3:
  488. #asm
  489.       FLDL2E
  490. #endasm
  491.        break;
  492.  
  493.      case 4:
  494. #asm
  495.       FLDL2T
  496. #endasm
  497.        break;
  498.  
  499.      case 5:
  500. #asm
  501.       FLDLG2
  502. #endasm
  503.        break;
  504.  
  505.      case 6:
  506. #asm
  507.       FLDLN2
  508. #endasm
  509.       break;
  510.     }
  511.     return;
  512. }
  513.  
  514. /*
  515. ** fltprint(n,num)
  516. ** Entry:
  517. **  num points to an 8-byte floating-point number
  518. **  n is the number of characters to print
  519. ** This function prints the floating-point number in scientific
  520. ** notation: -x.xxxxxx E-xx
  521. */
  522. fltprint(n,num) int n,num[];
  523. {
  524.     int ten[4];        /* Holder for ten */
  525.     char tmp[10];        /* Holder for one */
  526.     int msign;        /* Mantissa sign */
  527.     int powcnt;        /* Counter for power */
  528.     int i;            /* Loop count */
  529.     int j;            /* Counter */
  530.     int hld[4];        /* Holding val */
  531.  
  532. /*
  533. ** Copy the number we want to print out into a holding area,
  534. ** since we'll be altering it to print it out.
  535. */
  536.  
  537.     for(i=0;i<4;++i)
  538.         hld[i]=num[i];
  539.  
  540. /* First see what sign the number is.  Save in msign.
  541. ** If less than zero...change the sign of the number.
  542. ** If equal to zero, print it out and we can go home.
  543. */
  544.     msign=fcompz(hld);
  545.     if(msign==0) {
  546.       printf("0.0 E0");
  547.       return;
  548.     }
  549.     if(msign<0) {
  550.         fload(hld);
  551.         fabs();
  552.         fstore(hld);
  553.     }
  554.  
  555. /*
  556. ** Now see if the number is less than ten.  If so, multiply the number
  557. ** by ten repeatedly until it's not.  For each multiply, decrement a
  558. ** counter so we can keep track of the exponent.
  559. */
  560.     powcnt=0;
  561.     ten[0]=10;
  562.     ten[1]=ten[2]=ten[3]=0;
  563.     fint2float(ten);    /* Get floating-point 10 */
  564.     while(fcomp(hld,ten)==-1) {
  565.         --powcnt;
  566.         f2mult(hld,ten,hld);
  567.     }
  568.  
  569. /*
  570. ** Do the reverse of the above.  As long as the number is greater than or
  571. ** equal to 10, divide it by 10.  Increment the exponent counter for
  572. ** each multiply.
  573. */
  574.  
  575.     while(fcomp(hld,ten)!=-1) {
  576.         ++powcnt;
  577.         f2div(hld,ten,hld);
  578.     }
  579.  
  580. /*
  581. ** We now know the number is in the range 10 > n >= 1, and its
  582. ** exponent is in powcnt.  Print the number out as follows:
  583. ** Print out a minus sign if the number was less than 0.
  584. ** Multiply by 10 as needed to get the number of digits the
  585. ** caller wants printed to the left of the decimal point.
  586. ** Store the result in packed BCD format.
  587. ** Extract the BCD digits and print them out.
  588. */
  589.  
  590.     if(msign<0) printf("-");
  591.     fload(hld);
  592.     for(i=0;i<n-1;++i)
  593.         fmult(ten);
  594.     fstoreb(tmp);
  595.     for(i=0;i<n;++i) {
  596.         j=tmp[(n-i-1)/2];
  597.         if(((n-i)&1)==0) j=j>>4;
  598.         j=j&15;
  599.         printf("%d",j);
  600.         if(i==0) printf(".");
  601.     }
  602.  
  603. /* Now print out exponent */
  604.  
  605.     printf("E%d",powcnt);
  606.     return;
  607. }
  608.  
  609.