home *** CD-ROM | disk | FTP | other *** search
/ PC World 2000 August / PCWorld_2000-08_cd.bin / Software / TemaCD / xbasic / xbpro.exe / xb / xma.x < prev    next >
Text File  |  2000-01-12  |  57KB  |  1,965 lines
  1. '
  2. '
  3. ' ####################  Max Reason
  4. ' #####  PROLOG  #####  copyright 1988-2000
  5. ' ####################  XBasic mathematics function library
  6. '
  7. ' subject to GPLL license - see gpll.txt
  8. '
  9. ' maxreason@maxreason.com
  10. '
  11. ' for Windows XBasic
  12. ' for Linux XBasic
  13. '
  14. PROGRAM    "xma"
  15. VERSION    "0.0017"
  16. '
  17. IMPORT    "xst"
  18. '
  19. '
  20. EXPORT
  21. '
  22. ' ************************************
  23. ' *****  Math Library Functions  *****
  24. ' ************************************
  25. '
  26. ' Angles are always in RADIANS
  27. '
  28. DECLARE FUNCTION  Xma             ()
  29. DECLARE FUNCTION  XmaVersion$     ()
  30. DECLARE FUNCTION  DOUBLE  ACOS    (DOUBLE x)
  31. DECLARE FUNCTION  DOUBLE  ACOSH   (DOUBLE x)
  32. DECLARE FUNCTION  DOUBLE  ACOT    (DOUBLE x)
  33. DECLARE FUNCTION  DOUBLE  ACOTH   (DOUBLE x)
  34. DECLARE FUNCTION  DOUBLE  ACSC    (DOUBLE x)
  35. DECLARE FUNCTION  DOUBLE  ACSCH   (DOUBLE x)
  36. DECLARE FUNCTION  DOUBLE  ASEC    (DOUBLE x)
  37. DECLARE FUNCTION  DOUBLE  ASECH   (DOUBLE x)
  38. DECLARE FUNCTION  DOUBLE  ASIN    (DOUBLE x)
  39. DECLARE FUNCTION  DOUBLE  ASINH   (DOUBLE x)
  40. EXTERNAL FUNCTION  DOUBLE  ATAN    (DOUBLE x)
  41. DECLARE FUNCTION  DOUBLE  ATANH   (DOUBLE x)
  42. EXTERNAL FUNCTION  DOUBLE  COS     (DOUBLE x)
  43. DECLARE FUNCTION  DOUBLE  COSH    (DOUBLE x)
  44. DECLARE FUNCTION  DOUBLE  COT     (DOUBLE x)
  45. DECLARE FUNCTION  DOUBLE  COTH    (DOUBLE x)
  46. DECLARE FUNCTION  DOUBLE  CSC     (DOUBLE x)
  47. DECLARE FUNCTION  DOUBLE  CSCH    (DOUBLE x)
  48. EXTERNAL FUNCTION  DOUBLE  EXP     (DOUBLE x)
  49. DECLARE FUNCTION  DOUBLE  LOG     (DOUBLE x)
  50. EXTERNAL FUNCTION  DOUBLE  EXP2    (DOUBLE x)
  51. EXTERNAL FUNCTION  DOUBLE  EXP10   (DOUBLE x)
  52. DECLARE FUNCTION  DOUBLE  LOG10   (DOUBLE x)
  53. EXTERNAL FUNCTION  DOUBLE  POWER   (DOUBLE x, DOUBLE y)
  54. DECLARE FUNCTION  DOUBLE  SEC     (DOUBLE x)
  55. DECLARE FUNCTION  DOUBLE  SECH    (DOUBLE x)
  56. EXTERNAL FUNCTION  DOUBLE  SIN     (DOUBLE x)
  57. DECLARE FUNCTION  DOUBLE  SINH    (DOUBLE x)
  58. EXTERNAL FUNCTION  DOUBLE  SQRT    (DOUBLE x)
  59. EXTERNAL FUNCTION  DOUBLE  TAN     (DOUBLE x)
  60. DECLARE FUNCTION  DOUBLE  TANH    (DOUBLE x)
  61. '
  62. END EXPORT
  63. '
  64. INTERNAL FUNCTION  DOUBLE  Asin0   (DOUBLE x)
  65. INTERNAL FUNCTION  DOUBLE  Atan0   (DOUBLE x)
  66. INTERNAL FUNCTION  DOUBLE  Exp0    (DOUBLE x)
  67. INTERNAL FUNCTION  DOUBLE  Expmo   (DOUBLE x)
  68. INTERNAL FUNCTION  DOUBLE  Log0    (DOUBLE x)
  69. '
  70. ' For Xma:
  71. '
  72. INTERNAL FUNCTION  DOUBLE  Clog   (DOUBLE x)
  73. INTERNAL FUNCTION  DOUBLE  Clog0  (DOUBLE x)
  74. '
  75. EXTERNAL FUNCTION  XxxFCLEX ()
  76. EXTERNAL FUNCTION  XxxFINIT ()
  77. EXTERNAL FUNCTION  XxxFSTCW ()
  78. EXTERNAL FUNCTION  XxxFSTSW ()
  79. '
  80. EXTERNAL FUNCTION  DOUBLE  XxxF2XM1   (x#)
  81. EXTERNAL FUNCTION  DOUBLE  XxxFABS    (x#)
  82. EXTERNAL FUNCTION  DOUBLE  XxxFCHS    (x#)
  83. EXTERNAL FUNCTION  DOUBLE  XxxFCOS    (x#)
  84. EXTERNAL FUNCTION  DOUBLE  XxxFLDZ    ()
  85. EXTERNAL FUNCTION  DOUBLE  XxxFLD1    ()
  86. EXTERNAL FUNCTION  DOUBLE  XxxFLDPI   ()
  87. EXTERNAL FUNCTION  DOUBLE  XxxFLDL2E  ()
  88. EXTERNAL FUNCTION  DOUBLE  XxxFLDL2T  ()
  89. EXTERNAL FUNCTION  DOUBLE  XxxFLDLG2  ()
  90. EXTERNAL FUNCTION  DOUBLE  XxxFLDLN2  ()
  91. EXTERNAL FUNCTION  DOUBLE  XxxFPATAN  (x#, @y#)
  92. EXTERNAL FUNCTION  DOUBLE  XxxFPREM   (x#, @y#)
  93. EXTERNAL FUNCTION  DOUBLE  XxxFPREM1  (x#, @y#)
  94. EXTERNAL FUNCTION  DOUBLE  XxxFPTAN   (x#, @y#)
  95. EXTERNAL FUNCTION  DOUBLE  XxxFRNDINT (x#)
  96. EXTERNAL FUNCTION  DOUBLE  XxxFSCALE  (x#, y#)
  97. EXTERNAL FUNCTION  DOUBLE  XxxFSIN    (x#)
  98. EXTERNAL FUNCTION  DOUBLE  XxxFSINCOS (x#, @y#)
  99. EXTERNAL FUNCTION  DOUBLE  XxxFSQRT   (x#)
  100. EXTERNAL FUNCTION  DOUBLE  XxxFXTRACT (x#, @y#)
  101. EXTERNAL FUNCTION  DOUBLE  XxxFYL2X   (x#, @y#)
  102. EXTERNAL FUNCTION  DOUBLE  XxxFYL2XP1 (x#, @y#)
  103. '
  104. EXPORT
  105. '
  106. ' ************************************
  107. ' *****  Math Library Constants  *****
  108. ' ************************************
  109. '
  110.     $$NNAN            =    0dFFFFFFFFFFFFFFFF
  111.     $$PNAN            = 0d7FFFFFFFFFFFFFFF
  112.     $$NINF            = 0dFFF0000000000000
  113.     $$PINF            = 0d7FF0000000000000
  114.     $$RADIANS        = 1
  115.     $$DEGREES        = 2
  116.     $$DEGTORAD    = 0d3F91DF46A2529D39
  117.     $$RADTODEG    = 0d404CA5DC1A63C1F8
  118.     $$PI                = 0d400921FB54442D18
  119.     $$TWOPI            = 0d401921FB54442D18
  120.     $$PI3DIV2        = 0d4012D97C7F3321D2
  121.     $$PIDIV2        = 0d3FF921FB54442D18
  122.     $$PIDIV4        = 0d3FE921FB54442D18
  123.     $$INVPI            = 0d3FD45F306DC9C883
  124.     $$SQRT2            = 0d3FF6A09E667F3BCD
  125.     $$SQRT2DIV2    = 0d3FE6A09E667F3BCD
  126.     $$INVSQRT2    = 0d3FE6A09E667F3BCD
  127.     $$E                    = 0d4005BF0A8B145769
  128.     $$LOG2E            = 0d3FF71547652B82FE
  129.     $$LOG210        = 0d400A934F0979A371
  130.     $$LOGE2            = 0d3FE62E42FEFA39EF
  131.     $$LOGE10        = 0d4002EBB1BBB55516
  132.     $$LOGESQRT2    = 0d3FD62E42FEFA39EF
  133.     $$LOG102        = 0d3FD34413509F79FF
  134.     $$LOG10E        = 0d3FDBCB7B1526E50E
  135.     $$PIDIV8        = 0d3FD921FB54442D18
  136.     $$PI3DIV8        = 0d3FF2D97C7F3321D2
  137. END EXPORT
  138. '
  139. '
  140. ' ********************
  141. ' *****  Xma ()  *****
  142. ' ********************
  143. '
  144. FUNCTION  Xma ()
  145. '
  146.     a$ = "Max Reason"
  147.     a$ = "copyright 1988-2000"
  148.     a$ = "XBasic mathematics function library"
  149.     a$ = "maxreason@maxreason.com"
  150.     a$ = ""
  151. '
  152. '
  153. ' the following doesn't work unless the answer is less than 2
  154. ' because the stupid f2xm1 opcode only works if -1 <= x <= +1.
  155. '
  156. '    FOR x# = .01# TO 2# STEP .125#
  157. '        asw = XxxFSTSW ()
  158. '        acw = XxxFSTCW ()
  159. '        a# = EXP (x#)
  160. '        i# = XxxFETOX (x#)
  161. '        zsw = XxxFSTSW ()
  162. '        PRINT FORMAT$("###.###########",x#);;; FORMAT$("###.###########",a#);; FORMAT$("###.###########",i#);; HEX$(asw>>11,1);; HEX$(zsw>>11,1);; HEX$(acw,4)
  163. '    NEXT x#
  164. '    PRINT
  165. '
  166. ' the following doesn't work unless the answer is less than 2
  167. ' because the stupid f2xm1 opcode only works if -1 <= x <= +1.
  168. '
  169. '
  170. '    FOR x# = .01# TO 2# STEP .125#
  171. '        asw = XxxFSTSW ()
  172. '        b# = 10# ** x#
  173. '        j# = XxxFTENTOX (x#)
  174. '        zsw = XxxFSTSW ()
  175. '        PRINT FORMAT$("###.###########",x#);;; FORMAT$("###.###########",b#);; FORMAT$("###.###########",j#);; HEX$(asw>>11,1);; HEX$(zsw>>11,1)
  176. '    NEXT x#
  177. '    PRINT
  178. '
  179. ' the following doesn't work unless the answer is less than 2
  180. ' because the stupid f2xm1 opcode only works if -1 <= x <= +1.
  181. '
  182. '
  183. '    FOR x# = .1# TO 2.0001# STEP .02978#
  184. '        asw = XxxFSTSW ()
  185. '        c# = x# ** x#
  186. '        k# = XxxFYTOX (x#, x#)
  187. '        c$$ = GIANTAT (&c#)
  188. '        x$$ = GIANTAT (&x#)
  189. '        k$$ = GIANTAT (&k#)
  190. '        zsw = XxxFSTSW ()
  191. '        PRINT FORMAT$("###.###########",x#);;; FORMAT$("###.###########",c#);; FORMAT$("###.###########",k#);; HEX$(x$$,16);; HEX$(c$$,16);; HEX$(k$$,16);;; HEX$(asw>>11,1);; HEX$(zsw>>11,1)
  192. '    NEXT x#
  193. '    PRINT
  194.     RETURN
  195. '
  196. '
  197. ' test math library against pentium instructions
  198. '
  199. ' test constants
  200. '
  201. '    GOSUB TestFLDZ
  202. '    GOSUB TestFLD1
  203. '    GOSUB TestFLDPI
  204. '    GOSUB TestFLDL2E
  205. '    GOSUB TestFLDL2T
  206. '    GOSUB TestFLDLG2
  207. '    GOSUB TestFLDLN2
  208. '
  209. ' test functions
  210. '
  211. '    GOSUB TestF2XM1
  212. '    GOSUB TestFABS
  213. '    GOSUB TestFCHS
  214. '    GOSUB TestFCOS
  215. '    GOSUB TestFPATAN
  216. '    GOSUB TestFPREM
  217. '    GOSUB TestFPREM1
  218. '    GOSUB TestFPTAN
  219. '    GOSUB TestFRNDINT
  220. '    GOSUB TestFSCALE
  221. '    GOSUB TestFSIN
  222. '    GOSUB TestFSINCOS
  223. '    GOSUB TestFSQRT
  224. '    GOSUB TestFXTRACT
  225. '    GOSUB TestFYL2X
  226. '    GOSUB TestFYL2XP1
  227.     RETURN
  228. '
  229. '
  230. ' *****  test subroutines  *****
  231. '
  232. ' 0# vs XxxFLDZ
  233. '
  234. SUB TestFLDZ
  235.     a# = 0#
  236.     b# = XxxFLDZ ()
  237.     a$$ = GIANTAT(&a#)
  238.     b$$ = GIANTAT(&b#)
  239.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  240. END SUB
  241. '
  242. ' 1# vs XxxFLD1
  243. '
  244. SUB TestFLD1
  245.     a# = 1#
  246.     b# = XxxFLD1 ()
  247.     a$$ = GIANTAT(&a#)
  248.     b$$ = GIANTAT(&b#)
  249.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  250. END SUB
  251. '
  252. ' $$PI vs XxxFLDPI
  253. '
  254. SUB TestFLDPI
  255.     a# = $$PI
  256.     b# = XxxFLDPI ()
  257.     a$$ = GIANTAT(&a#)
  258.     b$$ = GIANTAT(&b#)
  259.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  260. END SUB
  261. '
  262. ' $$LOG2E vs XxxFLDL2E
  263. '
  264. SUB TestFLDL2E
  265.     a# = $$LOG2E
  266.     b# = XxxFLDL2E ()
  267.     a$$ = GIANTAT(&a#)
  268.     b$$ = GIANTAT(&b#)
  269.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  270. END SUB
  271. '
  272. ' $$LOG210 vs XxxFLDL2T
  273. '
  274. SUB TestFLDL2T
  275.     a# = $$LOG210
  276.     b# = XxxFLDL2T ()
  277.     a$$ = GIANTAT(&a#)
  278.     b$$ = GIANTAT(&b#)
  279.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  280. END SUB
  281. '
  282. ' $$LOG102 vs XxxFLDLG2
  283. '
  284. SUB TestFLDLG2
  285.     a# = $$LOG102
  286.     b# = XxxFLDLG2 ()
  287.     a$$ = GIANTAT(&a#)
  288.     b$$ = GIANTAT(&b#)
  289.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  290. END SUB
  291. '
  292. ' $$LOGE2 vs XxxFLDLN2
  293. '
  294. SUB TestFLDLN2
  295.     a# = $$LOGE2
  296.     b# = XxxFLDLN2 ()
  297.     a$$ = GIANTAT(&a#)
  298.     b$$ = GIANTAT(&b#)
  299.     PRINT HEX$(a$$,16);; HEX$(b$$,16)
  300. END SUB
  301. '
  302. ' (2# ** x# - 1#) vs XxxF2XM1
  303. '
  304. SUB TestF2XM1
  305.     FOR i# = -1# TO +1# STEP 1# / 64#
  306.         x# = 2# ** i# - 1#
  307.         y# = XxxF2XM1 (i#)
  308.         i$$ = GIANTAT (&i#)
  309.         x$$ = GIANTAT (&x#)
  310.         y$$ = GIANTAT (&y#)
  311.         PRINT FORMAT$ ("##.####", i#);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  312.     NEXT i#
  313. END SUB
  314. '
  315. ' ABS() vs XxxFABS
  316. '
  317. SUB TestFABS
  318.     FOR i# = -$$TWOPI TO $$TWOPI STEP $$TWOPI / 64#
  319.         x# = ABS (i#)
  320.         y# = XxxFABS (i#)
  321.         i$$ = GIANTAT (&i#)
  322.         x$$ = GIANTAT (&x#)
  323.         y$$ = GIANTAT (&y#)
  324.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  325.     NEXT i#
  326. END SUB
  327. '
  328. ' - vs XxxFCHS
  329. '
  330. SUB TestFCHS
  331.     FOR i# = -$$TWOPI TO $$TWOPI STEP $$TWOPI / 64#
  332.         x# = -i#
  333.         y# = XxxFCHS (i#)
  334.         i$$ = GIANTAT (&i#)
  335.         x$$ = GIANTAT (&x#)
  336.         y$$ = GIANTAT (&y#)
  337.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  338.     NEXT i#
  339. END SUB
  340. '
  341. ' COS() vs XxxFCOS
  342. '
  343. SUB TestFCOS
  344.     FOR i# = 0 TO $$TWOPI * 2# STEP $$TWOPI / 64#
  345.         x# = COS (i#)
  346.         y# = XxxFCOS (i#)
  347.         i$$ = GIANTAT (&i#)
  348.         x$$ = GIANTAT (&x#)
  349.         y$$ = GIANTAT (&y#)
  350.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  351.     NEXT i#
  352. END SUB
  353. '
  354. ' ATAN() vs XxxFPATAN
  355. '
  356. SUB TestFPATAN
  357.     FOR i# = 0# TO $$TWOPI * 2# STEP $$TWOPI / 64#
  358.         x# = ATAN (i#)
  359.         y# = XxxFPATAN (i#, 1#)
  360.         i$$ = GIANTAT (&i#)
  361.         x$$ = GIANTAT (&x#)
  362.         y$$ = GIANTAT (&y#)
  363.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  364.     NEXT i#
  365. END SUB
  366. '
  367. ' XxxFPREM
  368. '
  369. SUB TestFPREM
  370.     FOR i# = 10# TO 14# STEP .5#
  371.         FOR j# = -2# TO 2# STEP .35#
  372.             x# = XxxFPREM (i#, j#)
  373.             y# = XxxFPREM1 (i#, j#)
  374.             PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.####", j#);; FORMAT$ ("##.####", x#);; FORMAT$ ("##.####", y#)
  375.         NEXT j#
  376.     NEXT i#
  377. END SUB
  378. '
  379. ' XxxFPREM1
  380. '
  381. SUB TestFPREM1
  382.     FOR i# = 10# TO 14# STEP .5#
  383.         FOR j# = -2# TO 2# STEP .35#
  384.             x# = XxxFPREM1 (i#, j#)
  385.             y# = XxxFPREM (i#, j#)
  386.             PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.####", j#);; FORMAT$ ("##.####", x#);; FORMAT$ ("##.####", y#)
  387.         NEXT j#
  388.     NEXT i#
  389. END SUB
  390. '
  391. ' TAN() vs XxxFPTAN
  392. '
  393. SUB TestFPTAN
  394.     FOR i# = 0 TO $$TWOPI * 2# STEP $$TWOPI / 64#
  395.         x# = TAN (i#)
  396.         k# = XxxFPTAN (i#, @j#)
  397.         y# = k# / j#
  398.         i$$ = GIANTAT (&i#)
  399.         x$$ = GIANTAT (&x#)
  400.         y$$ = GIANTAT (&y#)
  401.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  402.     NEXT i#
  403. END SUB
  404. '
  405. ' INT() and FIX() vs XxxFRNDINT
  406. '
  407. SUB TestFRNDINT
  408.     FOR i# = -2# TO +2# STEP .25#
  409.         x# = INT (i#)
  410.         y# = FIX (i#)
  411.         z# = XxxFRNDINT (i#)
  412. '        i$$ = GIANTAT (&i#)
  413. '        x$$ = GIANTAT (&x#)
  414. '        y$$ = GIANTAT (&y#)
  415. '        z$$ = GIANTAT (&z#)
  416. '        PRINT FORMAT$ ("##.####", i#);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(z$$,16);; HEX$(z$$-y$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#);; FORMAT$ ("##.##################", z#)
  417.         PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#);; FORMAT$ ("##.##################", z#)
  418.     NEXT i#
  419. END SUB
  420. '
  421. ' (2# ** i * j#) vs XxxFSCALE
  422. '
  423. SUB TestFSCALE
  424.     FOR i# = -2# TO 2# STEP .25#
  425.         FOR j# = -2# TO 2# STEP .25#
  426.             x# = 2# ** DOUBLE(FIX(i#)) * j#
  427.             y# = XxxFSCALE (i#, j#)
  428.             i$$ = GIANTAT (&i#)
  429.             x$$ = GIANTAT (&x#)
  430.             y$$ = GIANTAT (&y#)
  431.             PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.####", j#);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  432.         NEXT j#
  433.     NEXT i#
  434. END SUB
  435. '
  436. ' SIN() vs XxxFSIN
  437. '
  438. SUB TestFSIN
  439.     FOR i# = 0 TO $$TWOPI * 2# STEP $$TWOPI / 64#
  440.         x# = SIN (i#)
  441.         y# = XxxFSIN (i#)
  442.         i$$ = GIANTAT (&i#)
  443.         x$$ = GIANTAT (&x#)
  444.         y$$ = GIANTAT (&y#)
  445.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  446.     NEXT i#
  447. END SUB
  448. '
  449. ' SIN() and COS() vs XxxFSINCOS
  450. '
  451. SUB TestFSINCOS
  452.     FOR i# = 0 TO $$TWOPI * 2# STEP $$TWOPI / 64#
  453.         x# = SIN (i#)
  454.         xx# = COS (i#)
  455.         y# = XxxFSINCOS (i#, @z#)
  456.         i$$ = GIANTAT (&i#)
  457.         x$$ = GIANTAT (&x#)
  458.         xx$$ = GIANTAT (&xx#)
  459.         y$$ = GIANTAT (&y#)
  460.         z$$ = GIANTAT (&z#)
  461.         PRINT FORMAT$ ("##.####", i#/$$PIDIV2);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(xx$$,16);; HEX$(y$$,16);; HEX$(z$$,16);; HEX$(y$$-x$$,16);; HEX$(z$$-xx$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  462.     NEXT i#
  463. END SUB
  464. '
  465. ' SQRT() vs XxxFSQRT
  466. '
  467. SUB TestFSQRT
  468.     FOR i# = 0 TO $$TWOPI * 2# STEP $$TWOPI / 64#
  469.         x# = SQRT (i#)
  470.         y# = XxxFSQRT (i#)
  471.         i$$ = GIANTAT (&i#)
  472.         x$$ = GIANTAT (&x#)
  473.         y$$ = GIANTAT (&y#)
  474.         PRINT FORMAT$ ("##.####", i#);; HEX$(i$$,16);; HEX$(x$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", x#);; FORMAT$ ("##.##################", y#)
  475.     NEXT i#
  476. END SUB
  477. '
  478. ' XxxFXTRACT
  479. '
  480. SUB TestFXTRACT
  481.     FOR i# = .25# TO 2# STEP 1# / 64#
  482.         y# = XxxFXTRACT (i#, @j#)
  483.         i$$ = GIANTAT (&i#)
  484.         j$$ = GIANTAT (&j#)
  485.         y$$ = GIANTAT (&y#)
  486.         PRINT FORMAT$ ("##.####", i#);; HEX$(i$$,16);; HEX$(j$$,16);; HEX$(y$$,16);; FORMAT$ ("##.##################", y#);; FORMAT$ ("##.##################", j#)
  487.     NEXT i#
  488. END SUB
  489. '
  490. ' XxxFYL2X
  491. '
  492. SUB TestFYL2X
  493.     FOR i# = 1# TO 2# STEP .25#
  494.         FOR j# = 1# TO 2# STEP .25#
  495.             y# = XxxFYL2X (i#, j#)
  496.             i$$ = GIANTAT (&i#)
  497.             j$$ = GIANTAT (&j#)
  498.             y$$ = GIANTAT (&y#)
  499.             z$$ = GIANTAT (&z#)
  500.             PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.####", j#);; HEX$(i$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", y#)
  501.         NEXT j#
  502.     NEXT i#
  503. END SUB
  504. '
  505. ' XxxFYL2XP1
  506. '
  507. SUB TestFYL2XP1
  508.     FOR i# = 1# TO 2# STEP .25#
  509.         FOR j# = 1# TO 2# STEP .25#
  510.             y# = XxxFYL2XP1 (i#, j#)
  511.             i$$ = GIANTAT (&i#)
  512.             j$$ = GIANTAT (&j#)
  513.             y$$ = GIANTAT (&y#)
  514.             z$$ = GIANTAT (&z#)
  515.             PRINT FORMAT$ ("##.####", i#);; FORMAT$ ("##.####", j#);; HEX$(i$$,16);; HEX$(y$$,16);; HEX$(y$$-x$$,16);; FORMAT$ ("##.##################", y#)
  516.         NEXT j#
  517.     NEXT i#
  518. END SUB
  519. END FUNCTION
  520. '
  521. '
  522. '  ****************************
  523. '  *****  XmaVersion$ ()  *****
  524. '  ****************************
  525. '
  526. FUNCTION  XmaVersion$ ()
  527.     version$ = VERSION$ (0)
  528.     RETURN (version$)
  529. END FUNCTION
  530. '
  531. '
  532. ' *********************
  533. ' *****  ACOS ()  *****
  534. ' *********************
  535. '
  536. FUNCTION DOUBLE ACOS (DOUBLE v) DOUBLE
  537.     XLONG    ##ERROR
  538. '
  539.     SELECT CASE TRUE
  540.         CASE (v < -1#)    : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  541.         CASE (v > 1#)        : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  542.         CASE (v = -1#)    : RETURN ($$PI)
  543.         CASE (v =  0#)    : RETURN ($$PIDIV2)
  544.         CASE (v =  1#)    : RETURN (0#)
  545.         CASE ELSE                : RETURN ($$PIDIV2 - ASIN(v))
  546.     END SELECT
  547. END FUNCTION
  548. '
  549. '
  550. ' **********************
  551. ' *****  ACOSH ()  *****
  552. ' **********************
  553. '
  554. FUNCTION DOUBLE ACOSH (DOUBLE v) DOUBLE
  555.     XLONG    ##ERROR
  556. '
  557.     IF (v < 1#) THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  558.     RETURN (LOG(v + SQRT(v*v - 1#)))
  559. END FUNCTION
  560. '
  561. '
  562. ' *********************
  563. ' *****  ACOT ()  *****
  564. ' *********************
  565. '
  566. FUNCTION DOUBLE ACOT (DOUBLE v) DOUBLE
  567.      RETURN ($$PIDIV2 - ATAN(v))
  568. END FUNCTION
  569. '
  570. '
  571. ' **********************
  572. ' *****  ACOTH ()  *****
  573. ' **********************
  574. '
  575. FUNCTION DOUBLE ACOTH (DOUBLE v) DOUBLE
  576.     XLONG    ##ERROR
  577. '
  578.     SELECT CASE TRUE
  579.         CASE (v < -1#)    : RETURN (ATANH(1#/v))
  580.         CASE (v > +1#)    : RETURN (ATANH(1#/v))
  581.         CASE ELSE                : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  582.     END SELECT
  583. END FUNCTION
  584. '
  585. '
  586. ' *********************
  587. ' *****  ACSC ()  *****
  588. ' *********************
  589. '
  590. FUNCTION DOUBLE ACSC (DOUBLE v) DOUBLE
  591.     XLONG    ##ERROR
  592. '
  593.     IF (ABS(v) < 1#) THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  594.     RETURN (ASIN(1#/v))
  595. END FUNCTION
  596. '
  597. '
  598. ' **********************
  599. ' *****  ACSCH ()  *****
  600. ' **********************
  601. '
  602. FUNCTION DOUBLE ACSCH (DOUBLE v) DOUBLE
  603.     XLONG    ##ERROR
  604. '
  605.     SELECT CASE TRUE
  606.         CASE (v = 0#)    : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  607.         CASE ELSE            : RETURN (ASINH(1#/v))
  608.     END SELECT
  609. END FUNCTION
  610. '
  611. '
  612. ' *********************
  613. ' *****  ASEC ()  *****
  614. ' *********************
  615. '
  616. FUNCTION DOUBLE ASEC (DOUBLE v) DOUBLE
  617.     XLONG    ##ERROR
  618. '
  619.     IF (ABS(v) < 1#) THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  620.     RETURN ($$PIDIV2 - ASIN(1#/v))
  621. END FUNCTION
  622. '
  623. '
  624. ' **********************
  625. ' *****  ASECH ()  *****
  626. ' **********************
  627. '
  628. FUNCTION DOUBLE ASECH (DOUBLE v) DOUBLE
  629.     XLONG    ##ERROR
  630. '
  631.     SELECT CASE TRUE
  632.         CASE (v <= 0#)    : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  633.         CASE (v > 1#)        : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  634.         CASE (v = 1)        : RETURN (0#)
  635.         CASE ELSE                : RETURN (ACOSH(1#/v))
  636.     END SELECT
  637. END FUNCTION
  638. '
  639. '
  640. ' *********************
  641. ' *****  ASIN ()  *****  Returns values between -PI/2 and +PI/2 inclusive
  642. ' *********************
  643. '
  644. FUNCTION DOUBLE ASIN (DOUBLE a) DOUBLE
  645.     XLONG    ##ERROR
  646. '
  647.     SELECT CASE TRUE
  648.         CASE (a < -1#)    : ##ERROR = $$ErrorNatureInvalidArgument: RETURN (0#)
  649.         CASE (a > 1#)        : ##ERROR = $$ErrorNatureInvalidArgument: RETURN (0#)
  650.         CASE (a = -1#)    : RETURN (-$$PIDIV2)
  651.         CASE (a = 0#)        : RETURN (0#)
  652.         CASE (a = 1#)        : RETURN ($$PIDIV2)
  653.         CASE (a > 0#)        : theSign = +1#
  654.         CASE ELSE                : theSign = -1# : a = -a
  655.     END SELECT
  656. '
  657.     IF (a <= $$INVSQRT2) THEN                                    ' a <= sin(45)
  658.         aa = a * a
  659.         x = a * Asin0 (aa)
  660.     ELSE
  661.         aa = (.5# - (a * .5#))
  662.         a = SQRT (aa)
  663.         x = $$PIDIV2 - (2# * a * Asin0 (aa))        ' x = 90 - 2 arcsin(a)
  664.     END IF
  665.     RETURN (x * theSign)
  666. END FUNCTION
  667. '
  668. '
  669. ' **********************
  670. ' *****  ASINH ()  *****
  671. ' **********************
  672. '
  673. FUNCTION DOUBLE ASINH (DOUBLE v) DOUBLE
  674.     IFZ v THEN RETURN (0#)
  675.     RETURN (LOG(v + SQRT(v*v + 1#)))
  676. END FUNCTION
  677. '
  678. '
  679. ' *********************
  680. ' *****  ATAN ()  *****  Hart #5034
  681. ' *********************
  682. '
  683. 'FUNCTION DOUBLE ATAN (DOUBLE v) DOUBLE
  684. '    STATIC first, w1, x1 , y1, w2, x2, y2, w3, x3, y3, w4, x4, y4
  685. '
  686. ' FPU routine
  687. '
  688. '    RETURN (XxxFPATAN(v,1#))
  689. '
  690. ' hand coded routine
  691. '
  692. '    IFZ first THEN GOSUB InitVars
  693. '
  694. ' Handle sign of value and result
  695. '
  696. '    SELECT CASE TRUE
  697. '        CASE (v > 0#)    : theSign = +1#
  698. '        CASE (v = 0#)    : RETURN (0#)
  699. '        CASE ELSE            : theSign = -1# : v = -v
  700. '    END SELECT
  701. '
  702. ' Different routines for different sub-quadrants
  703. '
  704. '    SELECT CASE TRUE
  705. '        CASE (v < .19891236737965800691#)        : RETURN (Atan0(v) * theSign)        ' <  tan(pi/16)
  706. '        CASE (v < .66817863791929892000#)        : w = w1: x = x1: y = y1                ' <  tan(3pi/16)
  707. '        CASE (v < .14966057626654890176d1)    : w = w2: x = x2: y = y2                ' <  tan(5pi/16)
  708. '        CASE (v < .50273394921258481045d1)    : w = w3: x = x3: y = y3                ' <  tan(7pi/16)
  709. '        CASE ELSE                                                        : w = w4:    x = x4: y = y4                ' <= tan(8pi/16)
  710. '    END SELECT
  711. '
  712. '    t = x - (y / (x + v))
  713. '    RETURN ((Atan0(t) + w) * theSign)
  714. '
  715. '
  716. ' *****  Initialize some variables  *****
  717. '
  718. 'SUB InitVars
  719. '    j#    = TAN ($$PIDIV8)                            ' *****  22.50 degrees  *****
  720. '    w1    = $$PIDIV8
  721. '    x1    = 1 / j#
  722. '    y1    = 1 + (1 / (j# * j#))
  723. '
  724. '    w2    = $$PIDIV4                                        ' *****  45.00 degrees  *****
  725. '    x2    = 1#
  726. '    y2    = 2#
  727. '
  728. '    j#    = TAN ($$PI3DIV8)                            ' *****  67.50 degrees  *****
  729. '    w3    = $$PI3DIV8
  730. '    x3    = 1 / j#
  731. '    y3    = 1 + (1 / (j# * j#))
  732. '
  733. '    w4     = $$PIDIV2                                        ' *****  90.00 degrees  *****
  734. '    x4    = 0
  735. '    y4    = 1
  736. '    first = $$TRUE
  737. 'END SUB
  738. 'END FUNCTION
  739. '
  740. '
  741. ' **********************
  742. ' *****  ATANH ()  *****
  743. ' **********************
  744. '
  745. FUNCTION DOUBLE ATANH (DOUBLE v) DOUBLE
  746.     XLONG    ##ERROR
  747. '
  748.     SELECT CASE TRUE
  749.         CASE (v = 0#)                : RETURN (0#)
  750.         CASE (ABS(v) < 1#)    : RETURN (LOG((1# + v) / (1# - v)) * 0.5#)
  751.         CASE ELSE                        : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  752.     END SELECT
  753. END FUNCTION
  754. '
  755. '
  756. ' ********************
  757. ' *****  COS ()  *****  Hart #3823
  758. ' ********************
  759. '
  760. 'FUNCTION DOUBLE COS (DOUBLE a) DOUBLE
  761. '
  762. ' FPU routine
  763. '
  764. '    RETURN (XxxFCOS(a))
  765. '
  766. ' hand coded routine
  767. '
  768. '    IF (a < -$$TWOPI) THEN                    ' more negative than -360 degrees
  769. '        r# = a / $$TWOPI                            ' r# = angle / 360 degrees
  770. '        i# = FIX (r#)                                    ' i# = integer part of r#
  771. '        j# = (r# - i#)                                ' j# = 0 to 1 (units of $$TWOPI)
  772. '        j# = j# + 1#                                    ' j# = ditto
  773. '        a = j# * $$TWOPI                            ' a  = 0 to -$$TWOPI
  774. '    END IF
  775. '    IF (a < 0) THEN a = a + $$TWOPI    ' a  = 0 to $$TWOPI = 0 to 360 degrees
  776. '
  777. '    IF (a > $$TWOPI) THEN                        ' more positive than +360 degrees
  778. '        r# = a / $$TWOPI                            ' r# = angle / 360 degrees
  779. '        i# = FIX (r#)                                    ' i# = integer part of r#
  780. '        a = (r# - i#)                                    ' a  = 0 to 1 (units of $$TWOPI)
  781. '        a = a * $$TWOPI                                ' a  = 0 to $$TWOPI
  782. '    END IF
  783. '
  784. '    fold into 0-90 with appropriate sign
  785. '
  786. '    SELECT CASE TRUE
  787. '        CASE (a <= $$PIDIV2)    : theSign = +1#
  788. '        CASE (a <= $$PI)            : theSign = -1#    : a = $$PI - a
  789. '        CASE (a <= $$PI3DIV2)    : theSign = -1#    : a = a - $$PI
  790. '        CASE ELSE                            : theSign = +1#    : a = $$TWOPI - a
  791. '    END SELECT
  792. '
  793. '    IF (a > $$PIDIV4) THEN
  794. '        x = SIN ($$PIDIV2 - a)                        ' use sin if a > 45
  795. '    ELSE
  796. '      a = a / $$PIDIV4
  797. '      aa = a * a
  798. '      x = aa * -.38577620372d-12 + .11500497024263d-9
  799. '        x = x * aa - .2461136382637005d-7
  800. '      x = x * aa + .359086044588581953d-5
  801. '      x = x * aa - .32599188692668755044d-3
  802. '      x = x * aa + .1585434424381541089754d-1
  803. '      x = x * aa - .30842513753404245242414
  804. '      x = x * aa + .99999999999999999996415
  805. '    END IF
  806. '    RETURN (x * theSign)
  807. 'END FUNCTION
  808. '
  809. '
  810. ' *********************
  811. ' *****  COSH ()  *****
  812. ' *********************
  813. '
  814. FUNCTION DOUBLE COSH (DOUBLE v) DOUBLE
  815. '
  816.     SELECT CASE TRUE
  817.         CASE (v = 0#)            : RETURN (1#)
  818.         CASE (v >= 19#)        : RETURN (EXP(v) * 0.5#)
  819.         CASE (v <= -19#)    : RETURN (EXP(-v) * 0.5#)
  820.         CASE ELSE                    : RETURN ((EXP(v) + EXP(-v)) * 0.5#)
  821.     END SELECT
  822. END FUNCTION
  823. '
  824. '
  825. ' ********************
  826. ' *****  COT ()  *****  Hart #4287 inverted
  827. ' ********************
  828. '
  829. FUNCTION DOUBLE COT (DOUBLE a) DOUBLE
  830. '
  831. ' FPU routine
  832. '
  833.     k# = XxxFPTAN (a, @j#)
  834.     RETURN (j# / k#)
  835. '
  836. ' hand coded routine
  837. '
  838. '    fold into 0 - 360
  839. '
  840.     DO WHILE (a < 0#)
  841.         a = a + $$TWOPI
  842.     LOOP
  843.     DO WHILE (a >= $$TWOPI)
  844.         a = a - $$TWOPI
  845.     LOOP
  846. '
  847. '    fold into 0-90 with appropriate sign
  848. '
  849.     SELECT CASE TRUE
  850.         CASE (a <= $$PIDIV2)    : theSign = +1#
  851.         CASE (a <= $$PI)            : theSign = -1# :    a = $$PI - a
  852.         CASE (a <= $$PI3DIV2)    : theSign = +1# :    a = a - $$PI
  853.         CASE ELSE                            : theSign = -1# :    a = $$TWOPI - a
  854.     END SELECT
  855. '
  856.     IF (a = 0#) THEN RETURN ($$PINF)
  857. '
  858.     a = a / $$PIDIV4
  859.     aa = a * a
  860.     x = aa * .1751083054422421995601107123d-1 - .279527948722905226792457575d2
  861.     x = x * aa + .6171941889092866899718078509d4
  862.     x = x * aa - .3498924461690337314553322577d6
  863.     x = x * aa + .413160917052212537541564775d7
  864.     y = aa - .4976002053470988434868766867d3                ' was aa * 1# - 497.xxxxx
  865.     y = y * aa + .5497880219885277215781502314d5
  866.     y = y * aa - .1527149650334576959585193088d7
  867.     y = y * aa + .5260528179299214050832459853d7
  868.     RETURN (y / (a * x)) * theSign
  869. END FUNCTION
  870. '
  871. '
  872. ' *********************
  873. ' *****  COTH ()  *****
  874. ' *********************
  875. '
  876. FUNCTION DOUBLE COTH (DOUBLE v) DOUBLE
  877.     XLONG    ##ERROR
  878. '
  879.     SELECT CASE TRUE
  880.         CASE (v = 0#)            : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  881.         CASE (v >= 19#)        : RETURN (1#)
  882.         CASE (v <= -19#)    : RETURN (-1#)
  883.         CASE ELSE                    : RETURN ((EXP(v) + EXP(-v)) / (EXP(v) - EXP(-v)))
  884.     END SELECT
  885. END FUNCTION
  886. '
  887. '
  888. ' ********************
  889. ' *****  CSC ()  *****  1 / SIN()
  890. ' ********************
  891. '
  892. FUNCTION DOUBLE CSC (DOUBLE a)
  893. '
  894.     x# = SIN(a)
  895.     IFZ x# THEN
  896.         RETURN ($$PINF)
  897.     ELSE
  898.         RETURN (1# / x#)
  899.     END IF
  900. END FUNCTION
  901. '
  902. '
  903. ' *********************
  904. ' *****  CSCH ()  *****
  905. ' *********************
  906. '
  907. FUNCTION DOUBLE CSCH (DOUBLE v) DOUBLE
  908.     XLONG    ##ERROR
  909. '
  910.     IFZ v THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  911.     RETURN (1# / SINH(v))
  912. END FUNCTION
  913. '
  914. '
  915. ' ********************
  916. ' *****  EXP ()  *****
  917. ' ********************
  918. '
  919. 'FUNCTION  DOUBLE  EXP (DOUBLE x) DOUBLE
  920. '    SLONG i, j, k, exp
  921. '    XLONG    ##ERROR
  922. '
  923. '    IFZ x THEN RETURN (1#)
  924. '
  925. '    IF (x < 0) THEN
  926. '        neg = $$TRUE
  927. '        x = -x
  928. '    ELSE
  929. '        neg = $$FALSE
  930. '    END IF
  931. '    y        = x * $$LOG2E                        ' y  = x * log2(e)
  932. '    i        = FIX (y)                                ' i  = INT (y)
  933. '    i#    = i                                            ' i# = INTpart of y
  934. '    z        = y - i#                                ' z  = FRACpart of y
  935. '    IF (z <= .5#) THEN
  936. '        r = Exp0 (z)
  937. '    ELSE
  938. '        r = Exp0 (z - .5#)
  939. '        r = r * $$SQRT2
  940. '    END IF
  941. '    j        = DHIGH (r)
  942. '    k        = DLOW (r)
  943. '    exp    = j {11, 20}
  944. '    exp = exp - 1023 + i
  945. '    exp = exp + 1023
  946. '
  947. '    IF (exp < 0) THEN
  948. '        ##ERROR = $$ErrorNatureOverflow
  949. '        IF neg THEN
  950. '            RETURN ($$PINF)
  951. '        ELSE
  952. '            RETURN (0#)
  953. '        END IF
  954. '    END IF
  955. '    IF (exp > 2047)
  956. '        ##ERROR = $$ErrorNatureOverflow
  957. '        IF neg THEN
  958. '            RETURN (0#)
  959. '        ELSE
  960. '            RETURN ($$PINF)
  961. '        END IF
  962. '    END IF
  963. '
  964. '    j = j AND 0x800FFFFF
  965. '    j = j OR (exp << 20)
  966. '    r = DMAKE (j, k)
  967. '    IF neg THEN
  968. '        RETURN (1/r)
  969. '    ELSE
  970. '        RETURN (r)
  971. '    END IF
  972. 'END FUNCTION
  973. '
  974. '
  975. ' ********************
  976. ' *****  LOG ()  *****
  977. ' ********************
  978. '
  979. FUNCTION DOUBLE LOG (DOUBLE v) DOUBLE
  980.     SLONG  upper,  lower,  exp,  exp0,  exp1,  exp2
  981.     XLONG  ##ERROR
  982. '
  983.     K1 = DMAKE (0xBF2BD010, 0x5C610CA9)
  984.     K2 = SMAKE (0x3F318000)
  985. '
  986.     SELECT CASE TRUE
  987.         CASE (v <= 0#)    : ##ERROR = $$ErrorNatureInvalidArgument: RETURN (0#)
  988.         CASE (v = 1#)        : RETURN (0#)
  989.     END SELECT
  990. '
  991.     upper    = DHIGH (v)
  992.     lower    = DLOW (v)
  993.     exp        = upper {11, 20}
  994.     exp0    = upper AND 0x800FFFFF
  995.     exp1    = exp0 OR 0x3FE00000
  996.     exp2    = exp - 1022
  997.     frac    = DMAKE (exp1, lower)                        ' frac is between 1/2 and 1
  998.     IF (frac >= $$INVSQRT2) THEN                    ' frac must be between 1/sqrt2 and sqrt2
  999.         z = (frac - 1#) / (frac + 1#)                ' it is!
  1000.     ELSE
  1001.         DEC exp2                                                        ' nope: dec the exponent, transfer to frac
  1002.         z = (frac - .5#) / (frac + .5#)            ' mult frac by 2 (clever...)
  1003.     END IF
  1004.     zp = z * Log0 (z)
  1005.     j = ((exp2 * K1) + zp) + (exp2 * K2)
  1006. '    j = exp2 * $$LOGE2 + zp
  1007.     RETURN (j)
  1008. END FUNCTION
  1009. '
  1010. '
  1011. ' **********************
  1012. ' *****  EXP10 ()  *****
  1013. ' **********************
  1014. '
  1015. 'FUNCTION DOUBLE EXP10 (DOUBLE v) DOUBLE
  1016. '    RETURN (POWER(10#,v))
  1017. 'END FUNCTION
  1018. '
  1019. '
  1020. ' **********************
  1021. ' *****  LOG10 ()  *****
  1022. ' **********************
  1023. '
  1024. FUNCTION DOUBLE LOG10 (DOUBLE v) DOUBLE
  1025.     RETURN (LOG (v) * $$LOG10E)
  1026. END FUNCTION
  1027. '
  1028. '
  1029. ' **********************
  1030. ' *****  POWER ()  *****
  1031. ' **********************
  1032. '
  1033. 'FUNCTION DOUBLE POWER (DOUBLE x, DOUBLE y) DOUBLE
  1034. '    XLONG  ##ERROR
  1035. '
  1036. '    SELECT CASE TRUE
  1037. '        CASE (y = 1#)                                            : RETURN (x)
  1038. '        CASE ((x = 0#) AND (y <= 0#))            : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  1039. '        CASE (x = 0#)                                            : RETURN (0#)
  1040. '        CASE (y = 0#)                                            : RETURN (1#)
  1041. '        CASE ((x < 0#) AND (FIX(y) != y))    : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  1042. '        CASE ELSE                                                    : RETURN (EXP(y * LOG(x)))
  1043. '    END SELECT
  1044. 'END FUNCTION
  1045. '
  1046. '
  1047. ' ********************
  1048. ' *****  SEC ()  *****  1 / COS()
  1049. ' ********************
  1050. '
  1051. FUNCTION DOUBLE SEC (DOUBLE a) DOUBLE
  1052. '
  1053.     x = COS(a)
  1054.     IFZ x THEN
  1055.         RETURN ($$PINF)
  1056.     ELSE
  1057.       RETURN (1# / x)
  1058.     END IF
  1059. END FUNCTION
  1060. '
  1061. '
  1062. ' *********************
  1063. ' *****  SECH ()  *****
  1064. ' *********************
  1065. '
  1066. FUNCTION DOUBLE SECH (DOUBLE v) DOUBLE
  1067.     IFZ v THEN RETURN (1#)
  1068.     RETURN (1# / COSH(v))
  1069. END FUNCTION
  1070. '
  1071. '
  1072. ' ********************
  1073. ' *****  SIN ()  *****  Hart #3043
  1074. ' ********************
  1075. '
  1076. 'FUNCTION DOUBLE SIN (DOUBLE a) DOUBLE
  1077. '
  1078. ' FPU routine
  1079. '
  1080. '    RETURN (XxxFSIN(a))
  1081. '
  1082. ' hand coded routine
  1083. '
  1084. '    IF (a < -$$TWOPI) THEN                    ' more negative than -360 degrees
  1085. '        r# = a / $$TWOPI                            ' r# = angle / 360 degrees
  1086. '        i# = FIX (r#)                                    ' i# = integer part of r#
  1087. '        j# = (r# - i#)                                ' j# = 0 to 1 (units of $$TWOPI)
  1088. '        j# = j# + 1#                                    ' j# = ditto
  1089. '        a = j# * $$TWOPI                            ' a  = 0 to -$$TWOPI
  1090. '    END IF
  1091. '    IF (a < 0) THEN a = a + $$TWOPI    ' a  = 0 to $$TWOPI = 0 to 360 degrees
  1092. '
  1093. '    IF (a > $$TWOPI) THEN                        ' more positive than +360 degrees
  1094. '        r# = a / $$TWOPI                            ' r# = angle / 360 degrees
  1095. '        i# = FIX (r#)                                    ' i# = integer part of r#
  1096. '        a = (r# - i#)                                    ' a  = 0 to 1 (units of $$TWOPI)
  1097. '        a = a * $$TWOPI                                ' a  = 0 to $$TWOPI
  1098. '    END IF
  1099. '
  1100. '    fold into 0 - 90 with appropriate sign
  1101. '
  1102. '    SELECT CASE TRUE
  1103. '        CASE (a <= $$PIDIV2):        theSign = +1#
  1104. '        CASE (a <= $$PI):                theSign = +1#:    a = $$PI - a
  1105. '        CASE (a <= $$PI3DIV2):    theSign = -1#:    a = a - $$PI
  1106. '        CASE ELSE:                            theSign = -1#:    a = $$TWOPI - a
  1107. '    END SELECT
  1108. '
  1109. '    IF (a > $$PIDIV4) THEN
  1110. '        x = COS ($$PIDIV2 - a)                                    ' use cos if a > 45
  1111. '    ELSE
  1112. '      a = a / $$PIDIV4
  1113. '      aa = a * a
  1114. '        x = aa * .6877100349d-11 - .1757149292755d-8
  1115. '        x = x * aa + .313361621661904d-6
  1116. '        x = x * aa - .36576204158455695d-4
  1117. '        x = x * aa + .2490394570188736117d-2
  1118. '        x = x * aa - .80745512188280530192d-1
  1119. '        x = x * aa + .785398163397448307014#
  1120. '        x = x * a
  1121. '    END IF
  1122. '    RETURN (x * theSign)
  1123. 'END FUNCTION
  1124. '
  1125. '
  1126. ' *********************
  1127. ' *****  SINH ()  *****  >>>  may need add'l stuff for exp(v) ~= exp(-v)
  1128. ' *********************
  1129. '
  1130. FUNCTION DOUBLE SINH (DOUBLE v) DOUBLE
  1131.     SELECT CASE TRUE
  1132.         CASE (v = 0#)            : RETURN (0#)
  1133.         CASE (v >= 19#)   : RETURN (EXP(v) * 0.5#)
  1134.         CASE (v <= -19#)    : RETURN (EXP(-v) * -0.5#)
  1135.         CASE (v > .3#)        : RETURN ((EXP(v) - EXP(-v)) * 0.5#)
  1136.         CASE (v < -.3#)        : RETURN ((EXP(v) - EXP(-v)) * 0.5#)
  1137.         CASE ELSE                    : dx = Expmo(v)
  1138.                                                 RETURN (0.5# * (dx + (dx / (dx + 1))))
  1139.     END SELECT
  1140. END FUNCTION
  1141. '
  1142. '
  1143. ' *********************
  1144. ' *****  SQRT ()  *****  Newton-Raphson Iteration
  1145. ' *********************
  1146. '
  1147. 'FUNCTION DOUBLE SQRT (DOUBLE v) DOUBLE
  1148. '    XLONG i, j, exp, exp0, exp1, exp2, xfrac, too
  1149. '    STATIC XLONG sqrt_table_odd[]
  1150. '    STATIC XLONG sqrt_table_even[]
  1151. '    XLONG  ##ERROR,  ##WHOMASK
  1152. '
  1153. '    IFZ v THEN RETURN (0#)
  1154. '    IF (v < 0#) THEN ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  1155. '
  1156. ' FPU routine
  1157. '
  1158. '    RETURN (XxxFSQRT(v))
  1159. '
  1160. ' hand coded routine
  1161. '
  1162. '    IFZ sqrt_table_odd[] THEN GOSUB FillEstimateArrays
  1163. '
  1164. '    exp0 = DHIGH(v)
  1165. '    exp1 = exp0{11, 20}
  1166. '    exp2 = exp1 - 1022
  1167. '    xfrac = exp0{8, 12}
  1168. '    IF (xfrac < 0) THEN PRINT "SQRT(): Error: (xfrac < 0)  "; xfrac : RETURN (0#)
  1169. '    IF (xfrac > 255) THEN PRINT "SQRT(): Error: (xfrac > 255)  "; xfrac : RETURN (0#)
  1170. '    IFZ exp2{1, 0} THEN
  1171. '        exp = (exp2 >>> 1) + 1022
  1172. '        exp = MAKE (exp, 20)
  1173. '        too = sqrt_table_even[xfrac]
  1174. '        exp = exp + too
  1175. '        e        = DMAKE (exp, 0)
  1176. '    ELSE
  1177. '        exp = (exp2 >>> 1) + 1023
  1178. '        exp = MAKE (exp, 20)
  1179. '        too = sqrt_table_odd[xfrac]
  1180. '        exp = exp + too
  1181. '        e        = DMAKE (exp, 0)
  1182. '    END IF
  1183. '    PRINT v; TAB (22); HEX$(DHIGH(v), 8);; HEX$(DLOW(v), 8); TAB (40); e; TAB (62); HEX$(DHIGH(e), 8);; HEX$(DLOW(e), 8)
  1184. '    FOR i = 0 TO 15
  1185. '        t = .5# * ((v / e) + e)
  1186. ''        PRINT e; TAB (22); HEX$(DHIGH(e), 8);; HEX$(DLOW(e), 8); TAB (40); t; TAB (62); HEX$(DHIGH(t), 8);; HEX$(DLOW(t), 8)
  1187. '        IF (e = t) THEN EXIT FOR
  1188. '        e = t
  1189. '    NEXT i
  1190. '    RETURN (e)
  1191. '
  1192. ' fill arrays that give starting estimate of sqrt(x) where exponent is odd/even
  1193. '
  1194. 'SUB FillEstimateArrays
  1195. '    hold = ##WHOMASK
  1196. '    ##WHOMASK = 0
  1197. '    DIM sqrt_table_odd[255]
  1198. '    DIM sqrt_table_even[255]
  1199. '    ##WHOMASK = hold
  1200. '    sqrt_table_odd[0x00]     = 0x00000000
  1201. '    sqrt_table_odd[0x01]     = 0x000007FE
  1202. '    sqrt_table_odd[0x02]     = 0x00000FF8
  1203. '    sqrt_table_odd[0x03]     = 0x000017EE
  1204. '    sqrt_table_odd[0x04]     = 0x00001FE0
  1205. '    sqrt_table_odd[0x05]     = 0x000027CE
  1206. '    sqrt_table_odd[0x06]     = 0x00002FB8
  1207. '    sqrt_table_odd[0x07]     = 0x0000379F
  1208. '    sqrt_table_odd[0x08]     = 0x00003F81
  1209. '    sqrt_table_odd[0x09]     = 0x00004760
  1210. '    sqrt_table_odd[0x0A]     = 0x00004F3B
  1211. '    sqrt_table_odd[0x0B]     = 0x00005713
  1212. '    sqrt_table_odd[0x0C]     = 0x00005EE6
  1213. '    sqrt_table_odd[0x0D]     = 0x000066B6
  1214. '    sqrt_table_odd[0x0E]     = 0x00006E82
  1215. '    sqrt_table_odd[0x0F]     = 0x0000764A
  1216. '    sqrt_table_odd[0x10]     = 0x00007E0F
  1217. '    sqrt_table_odd[0x11]     = 0x000085D0
  1218. '    sqrt_table_odd[0x12]     = 0x00008D8D
  1219. '    sqrt_table_odd[0x13]     = 0x00009547
  1220. '    sqrt_table_odd[0x14]     = 0x00009CFD
  1221. '    sqrt_table_odd[0x15]     = 0x0000A4B0
  1222. '    sqrt_table_odd[0x16]     = 0x0000AC5F
  1223. '    sqrt_table_odd[0x17]     = 0x0000B40B
  1224. '    sqrt_table_odd[0x18]     = 0x0000BBB3
  1225. '    sqrt_table_odd[0x19]     = 0x0000C357
  1226. '    sqrt_table_odd[0x1A]     = 0x0000CAF8
  1227. '    sqrt_table_odd[0x1B]     = 0x0000D296
  1228. '    sqrt_table_odd[0x1C]     = 0x0000DA30
  1229. '    sqrt_table_odd[0x1D]     = 0x0000E1C7
  1230. '    sqrt_table_odd[0x1E]     = 0x0000E95A
  1231. '    sqrt_table_odd[0x1F]     = 0x0000F0EA
  1232. '    sqrt_table_odd[0x20]     = 0x0000F876
  1233. '    sqrt_table_odd[0x21]     = 0x00010000
  1234. '    sqrt_table_odd[0x22]     = 0x00010785
  1235. '    sqrt_table_odd[0x23]     = 0x00010F08
  1236. '    sqrt_table_odd[0x24]     = 0x00011687
  1237. '    sqrt_table_odd[0x25]     = 0x00011E03
  1238. '    sqrt_table_odd[0x26]     = 0x0001257C
  1239. '    sqrt_table_odd[0x27]     = 0x00012CF1
  1240. '    sqrt_table_odd[0x28]     = 0x00013463
  1241. '    sqrt_table_odd[0x29]     = 0x00013BD2
  1242. '    sqrt_table_odd[0x2A]     = 0x0001433E
  1243. '    sqrt_table_odd[0x2B]     = 0x00014AA7
  1244. '    sqrt_table_odd[0x2C]     = 0x0001520C
  1245. '    sqrt_table_odd[0x2D]     = 0x0001596F
  1246. '    sqrt_table_odd[0x2E]     = 0x000160CE
  1247. '    sqrt_table_odd[0x2F]     = 0x0001682A
  1248. '    sqrt_table_odd[0x30]     = 0x00016F83
  1249. '    sqrt_table_odd[0x31]     = 0x000176D9
  1250. '    sqrt_table_odd[0x32]     = 0x00017E2B
  1251. '    sqrt_table_odd[0x33]     = 0x0001857B
  1252. '    sqrt_table_odd[0x34]     = 0x00018CC8
  1253. '    sqrt_table_odd[0x35]     = 0x00019411
  1254. '    sqrt_table_odd[0x36]     = 0x00019B58
  1255. '    sqrt_table_odd[0x37]     = 0x0001A29B
  1256. '    sqrt_table_odd[0x38]     = 0x0001A9DC
  1257. '    sqrt_table_odd[0x39]     = 0x0001B11A
  1258. '    sqrt_table_odd[0x3A]     = 0x0001B854
  1259. '    sqrt_table_odd[0x3B]     = 0x0001BF8C
  1260. '    sqrt_table_odd[0x3C]     = 0x0001C6C1
  1261. '    sqrt_table_odd[0x3D]     = 0x0001CDF3
  1262. '    sqrt_table_odd[0x3E]     = 0x0001D522
  1263. '    sqrt_table_odd[0x3F]     = 0x0001DC4E
  1264. '    sqrt_table_odd[0x40]     = 0x0001E377
  1265. '    sqrt_table_odd[0x41]     = 0x0001EA9D
  1266. '    sqrt_table_odd[0x42]     = 0x0001F1C1
  1267. '    sqrt_table_odd[0x43]     = 0x0001F8E2
  1268. '    sqrt_table_odd[0x44]     = 0x00020000
  1269. '    sqrt_table_odd[0x45]     = 0x0002071B
  1270. '    sqrt_table_odd[0x46]     = 0x00020E33
  1271. '    sqrt_table_odd[0x47]     = 0x00021548
  1272. '    sqrt_table_odd[0x48]     = 0x00021C5B
  1273. '    sqrt_table_odd[0x49]     = 0x0002236B
  1274. '    sqrt_table_odd[0x4A]     = 0x00022A78
  1275. '    sqrt_table_odd[0x4B]     = 0x00023183
  1276. '    sqrt_table_odd[0x4C]     = 0x0002388A
  1277. '    sqrt_table_odd[0x4D]     = 0x00023F8F
  1278. '    sqrt_table_odd[0x4E]     = 0x00024692
  1279. '    sqrt_table_odd[0x4F]     = 0x00024D91
  1280. '    sqrt_table_odd[0x50]     = 0x0002548E
  1281. '    sqrt_table_odd[0x51]     = 0x00025B89
  1282. '    sqrt_table_odd[0x52]     = 0x00026280
  1283. '    sqrt_table_odd[0x53]     = 0x00026975
  1284. '    sqrt_table_odd[0x54]     = 0x00027068
  1285. '    sqrt_table_odd[0x55]     = 0x00027757
  1286. '    sqrt_table_odd[0x56]     = 0x00027E45
  1287. '    sqrt_table_odd[0x57]     = 0x0002852F
  1288. '    sqrt_table_odd[0x58]     = 0x00028C17
  1289. '    sqrt_table_odd[0x59]     = 0x000292FD
  1290. '    sqrt_table_odd[0x5A]     = 0x000299E0
  1291. '    sqrt_table_odd[0x5B]     = 0x0002A0C0
  1292. '    sqrt_table_odd[0x5C]     = 0x0002A79E
  1293. '    sqrt_table_odd[0x5D]     = 0x0002AE79
  1294. '    sqrt_table_odd[0x5E]     = 0x0002B552
  1295. '    sqrt_table_odd[0x5F]     = 0x0002BC28
  1296. '    sqrt_table_odd[0x60]     = 0x0002C2FC
  1297. '    sqrt_table_odd[0x61]     = 0x0002C9CD
  1298. '    sqrt_table_odd[0x62]     = 0x0002D09C
  1299. '    sqrt_table_odd[0x63]     = 0x0002D768
  1300. '    sqrt_table_odd[0x64]     = 0x0002DE32
  1301. '    sqrt_table_odd[0x65]     = 0x0002E4FA
  1302. '    sqrt_table_odd[0x66]     = 0x0002EBBF
  1303. '    sqrt_table_odd[0x67]     = 0x0002F281
  1304. '    sqrt_table_odd[0x68]     = 0x0002F942
  1305. '    sqrt_table_odd[0x69]     = 0x00030000
  1306. '    sqrt_table_odd[0x6A]     = 0x000306BB
  1307. '    sqrt_table_odd[0x6B]     = 0x00030D74
  1308. '    sqrt_table_odd[0x6C]     = 0x0003142B
  1309. '    sqrt_table_odd[0x6D]     = 0x00031ADF
  1310. '    sqrt_table_odd[0x6E]     = 0x00032191
  1311. '    sqrt_table_odd[0x6F]     = 0x00032841
  1312. '    sqrt_table_odd[0x70]     = 0x00032EEE
  1313. '    sqrt_table_odd[0x71]     = 0x00033599
  1314. '    sqrt_table_odd[0x72]     = 0x00033C42
  1315. '    sqrt_table_odd[0x73]     = 0x000342E8
  1316. '    sqrt_table_odd[0x74]     = 0x0003498C
  1317. '    sqrt_table_odd[0x75]     = 0x0003502E
  1318. '    sqrt_table_odd[0x76]     = 0x000356CD
  1319. '    sqrt_table_odd[0x77]     = 0x00035D6B
  1320. '    sqrt_table_odd[0x78]     = 0x00036406
  1321. '    sqrt_table_odd[0x79]     = 0x00036A9E
  1322. '    sqrt_table_odd[0x7A]     = 0x00037135
  1323. '    sqrt_table_odd[0x7B]     = 0x000377C9
  1324. '    sqrt_table_odd[0x7C]     = 0x00037E5B
  1325. '    sqrt_table_odd[0x7D]     = 0x000384EB
  1326. '    sqrt_table_odd[0x7E]     = 0x00038B79
  1327. '    sqrt_table_odd[0x7F]     = 0x00039204
  1328. '    sqrt_table_odd[0x80]     = 0x0003988E
  1329. '    sqrt_table_odd[0x81]     = 0x00039F15
  1330. '    sqrt_table_odd[0x82]     = 0x0003A59A
  1331. '    sqrt_table_odd[0x83]     = 0x0003AC1C
  1332. '    sqrt_table_odd[0x84]     = 0x0003B29D
  1333. '    sqrt_table_odd[0x85]     = 0x0003B91B
  1334. '    sqrt_table_odd[0x86]     = 0x0003BF98
  1335. '    sqrt_table_odd[0x87]     = 0x0003C612
  1336. '    sqrt_table_odd[0x88]     = 0x0003CC8A
  1337. '    sqrt_table_odd[0x89]     = 0x0003D300
  1338. '    sqrt_table_odd[0x8A]     = 0x0003D974
  1339. '    sqrt_table_odd[0x8B]     = 0x0003DFE6
  1340. '    sqrt_table_odd[0x8C]     = 0x0003E655
  1341. '    sqrt_table_odd[0x8D]     = 0x0003ECC3
  1342. '    sqrt_table_odd[0x8E]     = 0x0003F32F
  1343. '    sqrt_table_odd[0x8F]     = 0x0003F998
  1344. '    sqrt_table_odd[0x90]     = 0x00040000
  1345. '    sqrt_table_odd[0x91]     = 0x00040665
  1346. '    sqrt_table_odd[0x92]     = 0x00040CC8
  1347. '    sqrt_table_odd[0x93]     = 0x0004132A
  1348. '    sqrt_table_odd[0x94]     = 0x00041989
  1349. '    sqrt_table_odd[0x95]     = 0x00041FE6
  1350. '    sqrt_table_odd[0x96]     = 0x00042641
  1351. '    sqrt_table_odd[0x97]     = 0x00042C9B
  1352. '    sqrt_table_odd[0x98]     = 0x000432F2
  1353. '    sqrt_table_odd[0x99]     = 0x00043947
  1354. '    sqrt_table_odd[0x9A]     = 0x00043F9A
  1355. '    sqrt_table_odd[0x9B]     = 0x000445EC
  1356. '    sqrt_table_odd[0x9C]     = 0x00044C3B
  1357. '    sqrt_table_odd[0x9D]     = 0x00045288
  1358. '    sqrt_table_odd[0x9E]     = 0x000458D4
  1359. '    sqrt_table_odd[0x9F]     = 0x00045F1D
  1360. '    sqrt_table_odd[0xA0]     = 0x00046565
  1361. '    sqrt_table_odd[0xA1]     = 0x00046BAA
  1362. '    sqrt_table_odd[0xA2]     = 0x000471EE
  1363. '    sqrt_table_odd[0xA3]     = 0x00047830
  1364. '    sqrt_table_odd[0xA4]     = 0x00047E70
  1365. '    sqrt_table_odd[0xA5]     = 0x000484AE
  1366. '    sqrt_table_odd[0xA6]     = 0x00048AEA
  1367. '    sqrt_table_odd[0xA7]     = 0x00049124
  1368. '    sqrt_table_odd[0xA8]     = 0x0004975C
  1369. '    sqrt_table_odd[0xA9]     = 0x00049D93
  1370. '    sqrt_table_odd[0xAA]     = 0x0004A3C7
  1371. '    sqrt_table_odd[0xAB]     = 0x0004A9FA
  1372. '    sqrt_table_odd[0xAC]     = 0x0004B02B
  1373. '    sqrt_table_odd[0xAD]     = 0x0004B65A
  1374. '    sqrt_table_odd[0xAE]     = 0x0004BC87
  1375. '    sqrt_table_odd[0xAF]     = 0x0004C2B2
  1376. '    sqrt_table_odd[0xB0]     = 0x0004C8DC
  1377. '    sqrt_table_odd[0xB1]     = 0x0004CF03
  1378. '    sqrt_table_odd[0xB2]     = 0x0004D529
  1379. '    sqrt_table_odd[0xB3]     = 0x0004DB4D
  1380. '    sqrt_table_odd[0xB4]     = 0x0004E16F
  1381. '    sqrt_table_odd[0xB5]     = 0x0004E790
  1382. '    sqrt_table_odd[0xB6]     = 0x0004EDAE
  1383. '    sqrt_table_odd[0xB7]     = 0x0004F3CB
  1384. '    sqrt_table_odd[0xB8]     = 0x0004F9E6
  1385. '    sqrt_table_odd[0xB9]     = 0x00050000
  1386. '    sqrt_table_odd[0xBA]     = 0x00050617
  1387. '    sqrt_table_odd[0xBB]     = 0x00050C2D
  1388. '    sqrt_table_odd[0xBC]     = 0x00051241
  1389. '    sqrt_table_odd[0xBD]     = 0x00051853
  1390. '    sqrt_table_odd[0xBE]     = 0x00051E63
  1391. '    sqrt_table_odd[0xBF]     = 0x00052472
  1392. '    sqrt_table_odd[0xC0]     = 0x00052A7F
  1393. '    sqrt_table_odd[0xC1]     = 0x0005308A
  1394. '    sqrt_table_odd[0xC2]     = 0x00053694
  1395. '    sqrt_table_odd[0xC3]     = 0x00053C9C
  1396. '    sqrt_table_odd[0xC4]     = 0x000542A2
  1397. '    sqrt_table_odd[0xC5]     = 0x000548A6
  1398. '    sqrt_table_odd[0xC6]     = 0x00054EA9
  1399. '    sqrt_table_odd[0xC7]     = 0x000554AA
  1400. '    sqrt_table_odd[0xC8]     = 0x00055AAA
  1401. '    sqrt_table_odd[0xC9]     = 0x000560A7
  1402. '    sqrt_table_odd[0xCA]     = 0x000566A3
  1403. '    sqrt_table_odd[0xCB]     = 0x00056C9D
  1404. '    sqrt_table_odd[0xCC]     = 0x00057296
  1405. '    sqrt_table_odd[0xCD]     = 0x0005788D
  1406. '    sqrt_table_odd[0xCE]     = 0x00057E82
  1407. '    sqrt_table_odd[0xCF]     = 0x00058476
  1408. '    sqrt_table_odd[0xD0]     = 0x00058A68
  1409. '    sqrt_table_odd[0xD1]     = 0x00059059
  1410. '    sqrt_table_odd[0xD2]     = 0x00059647
  1411. '    sqrt_table_odd[0xD3]     = 0x00059C34
  1412. '    sqrt_table_odd[0xD4]     = 0x0005A220
  1413. '    sqrt_table_odd[0xD5]     = 0x0005A80A
  1414. '    sqrt_table_odd[0xD6]     = 0x0005ADF2
  1415. '    sqrt_table_odd[0xD7]     = 0x0005B3D9
  1416. '    sqrt_table_odd[0xD8]     = 0x0005B9BE
  1417. '    sqrt_table_odd[0xD9]     = 0x0005BFA1
  1418. '    sqrt_table_odd[0xDA]     = 0x0005C583
  1419. '    sqrt_table_odd[0xDB]     = 0x0005CB64
  1420. '    sqrt_table_odd[0xDC]     = 0x0005D142
  1421. '    sqrt_table_odd[0xDD]     = 0x0005D71F
  1422. '    sqrt_table_odd[0xDE]     = 0x0005DCFB
  1423. '    sqrt_table_odd[0xDF]     = 0x0005E2D5
  1424. '    sqrt_table_odd[0xE0]     = 0x0005E8AD
  1425. '    sqrt_table_odd[0xE1]     = 0x0005EE84
  1426. '    sqrt_table_odd[0xE2]     = 0x0005F45A
  1427. '    sqrt_table_odd[0xE3]     = 0x0005FA2D
  1428. '    sqrt_table_odd[0xE4]     = 0x00060000
  1429. '    sqrt_table_odd[0xE5]     = 0x000605D0
  1430. '    sqrt_table_odd[0xE6]     = 0x00060B9F
  1431. '    sqrt_table_odd[0xE7]     = 0x0006116D
  1432. '    sqrt_table_odd[0xE8]     = 0x00061739
  1433. '    sqrt_table_odd[0xE9]     = 0x00061D04
  1434. '    sqrt_table_odd[0xEA]     = 0x000622CD
  1435. '    sqrt_table_odd[0xEB]     = 0x00062894
  1436. '    sqrt_table_odd[0xEC]     = 0x00062E5A
  1437. '    sqrt_table_odd[0xED]     = 0x0006341F
  1438. '    sqrt_table_odd[0xEE]     = 0x000639E2
  1439. '    sqrt_table_odd[0xEF]     = 0x00063FA3
  1440. '    sqrt_table_odd[0xF0]     = 0x00064564
  1441. '    sqrt_table_odd[0xF1]     = 0x00064B22
  1442. '    sqrt_table_odd[0xF2]     = 0x000650DF
  1443. '    sqrt_table_odd[0xF3]     = 0x0006569B
  1444. '    sqrt_table_odd[0xF4]     = 0x00065C55
  1445. '    sqrt_table_odd[0xF5]     = 0x0006620E
  1446. '    sqrt_table_odd[0xF6]     = 0x000667C5
  1447. '    sqrt_table_odd[0xF7]     = 0x00066D7B
  1448. '    sqrt_table_odd[0xF8]     = 0x0006732F
  1449. '    sqrt_table_odd[0xF9]     = 0x000678E2
  1450. '    sqrt_table_odd[0xFA]     = 0x00067E93
  1451. '    sqrt_table_odd[0xFB]     = 0x00068443
  1452. '    sqrt_table_odd[0xFC]     = 0x000689F2
  1453. '    sqrt_table_odd[0xFD]     = 0x00068F9F
  1454. '    sqrt_table_odd[0xFE]     = 0x0006954B
  1455. '    sqrt_table_odd[0xFF]     = 0x00069AF5
  1456. '    sqrt_table_even[0x00]     = 0x0006A09E
  1457. '    sqrt_table_even[0x01]     = 0x0006ABEB
  1458. '    sqrt_table_even[0x02]     = 0x0006B733
  1459. '    sqrt_table_even[0x03]     = 0x0006C276
  1460. '    sqrt_table_even[0x04]     = 0x0006CDB2
  1461. '    sqrt_table_even[0x05]     = 0x0006D8E9
  1462. '    sqrt_table_even[0x06]     = 0x0006E41B
  1463. '    sqrt_table_even[0x07]     = 0x0006EF47
  1464. '    sqrt_table_even[0x08]     = 0x0006FA6E
  1465. '    sqrt_table_even[0x09]     = 0x00070590
  1466. '    sqrt_table_even[0x0A]     = 0x000710AC
  1467. '    sqrt_table_even[0x0B]     = 0x00071BC2
  1468. '    sqrt_table_even[0x0C]     = 0x000726D4
  1469. '    sqrt_table_even[0x0D]     = 0x000731E0
  1470. '    sqrt_table_even[0x0E]     = 0x00073CE7
  1471. '    sqrt_table_even[0x0F]     = 0x000747E8
  1472. '    sqrt_table_even[0x10]     = 0x000752E5
  1473. '    sqrt_table_even[0x11]     = 0x00075DDC
  1474. '    sqrt_table_even[0x12]     = 0x000768CE
  1475. '    sqrt_table_even[0x13]     = 0x000773BB
  1476. '    sqrt_table_even[0x14]     = 0x00077EA3
  1477. '    sqrt_table_even[0x15]     = 0x00078986
  1478. '    sqrt_table_even[0x16]     = 0x00079464
  1479. '    sqrt_table_even[0x17]     = 0x00079F3C
  1480. '    sqrt_table_even[0x18]     = 0x0007AA10
  1481. '    sqrt_table_even[0x19]     = 0x0007B4DF
  1482. '    sqrt_table_even[0x1A]     = 0x0007BFA9
  1483. '    sqrt_table_even[0x1B]     = 0x0007CA6E
  1484. '    sqrt_table_even[0x1C]     = 0x0007D52F
  1485. '    sqrt_table_even[0x1D]     = 0x0007DFEA
  1486. '    sqrt_table_even[0x1E]     = 0x0007EAA1
  1487. '    sqrt_table_even[0x1F]     = 0x0007F552
  1488. '    sqrt_table_even[0x20]     = 0x00080000
  1489. '    sqrt_table_even[0x21]     = 0x00080AA8
  1490. '    sqrt_table_even[0x22]     = 0x0008154B
  1491. '    sqrt_table_even[0x23]     = 0x00081FEA
  1492. '    sqrt_table_even[0x24]     = 0x00082A85
  1493. '    sqrt_table_even[0x25]     = 0x0008351A
  1494. '    sqrt_table_even[0x26]     = 0x00083FAB
  1495. '    sqrt_table_even[0x27]     = 0x00084A37
  1496. '    sqrt_table_even[0x28]     = 0x000854BF
  1497. '    sqrt_table_even[0x29]     = 0x00085F42
  1498. '    sqrt_table_even[0x2A]     = 0x000869C1
  1499. '    sqrt_table_even[0x2B]     = 0x0008743B
  1500. '    sqrt_table_even[0x2C]     = 0x00087EB1
  1501. '    sqrt_table_even[0x2D]     = 0x00088922
  1502. '    sqrt_table_even[0x2E]     = 0x0008938F
  1503. '    sqrt_table_even[0x2F]     = 0x00089DF8
  1504. '    sqrt_table_even[0x30]     = 0x0008A85C
  1505. '    sqrt_table_even[0x31]     = 0x0008B2BB
  1506. '    sqrt_table_even[0x32]     = 0x0008BD17
  1507. '    sqrt_table_even[0x33]     = 0x0008C76E
  1508. '    sqrt_table_even[0x34]     = 0x0008D1C0
  1509. '    sqrt_table_even[0x35]     = 0x0008DC0F
  1510. '    sqrt_table_even[0x36]     = 0x0008E659
  1511. '    sqrt_table_even[0x37]     = 0x0008F09F
  1512. '    sqrt_table_even[0x38]     = 0x0008FAE0
  1513. '    sqrt_table_even[0x39]     = 0x0009051E
  1514. '    sqrt_table_even[0x3A]     = 0x00090F57
  1515. '    sqrt_table_even[0x3B]     = 0x0009198C
  1516. '    sqrt_table_even[0x3C]     = 0x000923BD
  1517. '    sqrt_table_even[0x3D]     = 0x00092DEA
  1518. '    sqrt_table_even[0x3E]     = 0x00093813
  1519. '    sqrt_table_even[0x3F]     = 0x00094237
  1520. '    sqrt_table_even[0x40]     = 0x00094C58
  1521. '    sqrt_table_even[0x41]     = 0x00095674
  1522. '    sqrt_table_even[0x42]     = 0x0009608D
  1523. '    sqrt_table_even[0x43]     = 0x00096AA1
  1524. '    sqrt_table_even[0x44]     = 0x000974B2
  1525. '    sqrt_table_even[0x45]     = 0x00097EBE
  1526. '    sqrt_table_even[0x46]     = 0x000988C7
  1527. '    sqrt_table_even[0x47]     = 0x000992CB
  1528. '    sqrt_table_even[0x48]     = 0x00099CCC
  1529. '    sqrt_table_even[0x49]     = 0x0009A6C9
  1530. '    sqrt_table_even[0x4A]     = 0x0009B0C2
  1531. '    sqrt_table_even[0x4B]     = 0x0009BAB7
  1532. '    sqrt_table_even[0x4C]     = 0x0009C4A8
  1533. '    sqrt_table_even[0x4D]     = 0x0009CE95
  1534. '    sqrt_table_even[0x4E]     = 0x0009D87F
  1535. '    sqrt_table_even[0x4F]     = 0x0009E265
  1536. '    sqrt_table_even[0x50]     = 0x0009EC47
  1537. '    sqrt_table_even[0x51]     = 0x0009F625
  1538. '    sqrt_table_even[0x52]     = 0x000A0000
  1539. '    sqrt_table_even[0x53]     = 0x000A09D6
  1540. '    sqrt_table_even[0x54]     = 0x000A13A9
  1541. '    sqrt_table_even[0x55]     = 0x000A1D79
  1542. '    sqrt_table_even[0x56]     = 0x000A2744
  1543. '    sqrt_table_even[0x57]     = 0x000A310C
  1544. '    sqrt_table_even[0x58]     = 0x000A3AD1
  1545. '    sqrt_table_even[0x59]     = 0x000A4491
  1546. '    sqrt_table_even[0x5A]     = 0x000A4E4E
  1547. '    sqrt_table_even[0x5B]     = 0x000A5808
  1548. '    sqrt_table_even[0x5C]     = 0x000A61BE
  1549. '    sqrt_table_even[0x5D]     = 0x000A6B70
  1550. '    sqrt_table_even[0x5E]     = 0x000A751F
  1551. '    sqrt_table_even[0x5F]     = 0x000A7ECA
  1552. '    sqrt_table_even[0x60]     = 0x000A8872
  1553. '    sqrt_table_even[0x61]     = 0x000A9216
  1554. '    sqrt_table_even[0x62]     = 0x000A9BB7
  1555. '    sqrt_table_even[0x63]     = 0x000AA554
  1556. '    sqrt_table_even[0x64]     = 0x000AAEEE
  1557. '    sqrt_table_even[0x65]     = 0x000AB884
  1558. '    sqrt_table_even[0x66]     = 0x000AC217
  1559. '    sqrt_table_even[0x67]     = 0x000ACBA7
  1560. '    sqrt_table_even[0x68]     = 0x000AD533
  1561. '    sqrt_table_even[0x69]     = 0x000ADEBC
  1562. '    sqrt_table_even[0x6A]     = 0x000AE841
  1563. '    sqrt_table_even[0x6B]     = 0x000AF1C3
  1564. '    sqrt_table_even[0x6C]     = 0x000AFB41
  1565. '    sqrt_table_even[0x6D]     = 0x000B04BD
  1566. '    sqrt_table_even[0x6E]     = 0x000B0E35
  1567. '    sqrt_table_even[0x6F]     = 0x000B17A9
  1568. '    sqrt_table_even[0x70]     = 0x000B211B
  1569. '    sqrt_table_even[0x71]     = 0x000B2A89
  1570. '    sqrt_table_even[0x72]     = 0x000B33F3
  1571. '    sqrt_table_even[0x73]     = 0x000B3D5B
  1572. '    sqrt_table_even[0x74]     = 0x000B46BF
  1573. '    sqrt_table_even[0x75]     = 0x000B5020
  1574. '    sqrt_table_even[0x76]     = 0x000B597E
  1575. '    sqrt_table_even[0x77]     = 0x000B62D9
  1576. '    sqrt_table_even[0x78]     = 0x000B6C30
  1577. '    sqrt_table_even[0x79]     = 0x000B7584
  1578. '    sqrt_table_even[0x7A]     = 0x000B7ED6
  1579. '    sqrt_table_even[0x7B]     = 0x000B8824
  1580. '    sqrt_table_even[0x7C]     = 0x000B916E
  1581. '    sqrt_table_even[0x7D]     = 0x000B9AB6
  1582. '    sqrt_table_even[0x7E]     = 0x000BA3FB
  1583. '    sqrt_table_even[0x7F]     = 0x000BAD3C
  1584. '    sqrt_table_even[0x80]     = 0x000BB67A
  1585. '    sqrt_table_even[0x81]     = 0x000BBFB6
  1586. '    sqrt_table_even[0x82]     = 0x000BC8EE
  1587. '    sqrt_table_even[0x83]     = 0x000BD223
  1588. '    sqrt_table_even[0x84]     = 0x000BDB55
  1589. '    sqrt_table_even[0x85]     = 0x000BE484
  1590. '    sqrt_table_even[0x86]     = 0x000BEDB0
  1591. '    sqrt_table_even[0x87]     = 0x000BF6D9
  1592. '    sqrt_table_even[0x88]     = 0x000C0000
  1593. '    sqrt_table_even[0x89]     = 0x000C0923
  1594. '    sqrt_table_even[0x8A]     = 0x000C1243
  1595. '    sqrt_table_even[0x8B]     = 0x000C1B60
  1596. '    sqrt_table_even[0x8C]     = 0x000C247A
  1597. '    sqrt_table_even[0x8D]     = 0x000C2D91
  1598. '    sqrt_table_even[0x8E]     = 0x000C36A6
  1599. '    sqrt_table_even[0x8F]     = 0x000C3FB7
  1600. '    sqrt_table_even[0x90]     = 0x000C48C6
  1601. '    sqrt_table_even[0x91]     = 0x000C51D1
  1602. '    sqrt_table_even[0x92]     = 0x000C5ADA
  1603. '    sqrt_table_even[0x93]     = 0x000C63E0
  1604. '    sqrt_table_even[0x94]     = 0x000C6CE3
  1605. '    sqrt_table_even[0x95]     = 0x000C75E3
  1606. '    sqrt_table_even[0x96]     = 0x000C7EE0
  1607. '    sqrt_table_even[0x97]     = 0x000C87DA
  1608. '    sqrt_table_even[0x98]     = 0x000C90D2
  1609. '    sqrt_table_even[0x99]     = 0x000C99C7
  1610. '    sqrt_table_even[0x9A]     = 0x000CA2B9
  1611. '    sqrt_table_even[0x9B]     = 0x000CABA8
  1612. '    sqrt_table_even[0x9C]     = 0x000CB495
  1613. '    sqrt_table_even[0x9D]     = 0x000CBD7E
  1614. '    sqrt_table_even[0x9E]     = 0x000CC665
  1615. '    sqrt_table_even[0x9F]     = 0x000CCF49
  1616. '    sqrt_table_even[0xA0]     = 0x000CD82B
  1617. '    sqrt_table_even[0xA1]     = 0x000CE109
  1618. '    sqrt_table_even[0xA2]     = 0x000CE9E5
  1619. '    sqrt_table_even[0xA3]     = 0x000CF2BF
  1620. '    sqrt_table_even[0xA4]     = 0x000CFB95
  1621. '    sqrt_table_even[0xA5]     = 0x000D0469
  1622. '    sqrt_table_even[0xA6]     = 0x000D0D3A
  1623. '    sqrt_table_even[0xA7]     = 0x000D1609
  1624. '    sqrt_table_even[0xA8]     = 0x000D1ED5
  1625. '    sqrt_table_even[0xA9]     = 0x000D279E
  1626. '    sqrt_table_even[0xAA]     = 0x000D3064
  1627. '    sqrt_table_even[0xAB]     = 0x000D3928
  1628. '    sqrt_table_even[0xAC]     = 0x000D41EA
  1629. '    sqrt_table_even[0xAD]     = 0x000D4AA8
  1630. '    sqrt_table_even[0xAE]     = 0x000D5364
  1631. '    sqrt_table_even[0xAF]     = 0x000D5C1E
  1632. '    sqrt_table_even[0xB0]     = 0x000D64D5
  1633. '    sqrt_table_even[0xB1]     = 0x000D6D89
  1634. '    sqrt_table_even[0xB2]     = 0x000D763B
  1635. '    sqrt_table_even[0xB3]     = 0x000D7EEA
  1636. '    sqrt_table_even[0xB4]     = 0x000D8796
  1637. '    sqrt_table_even[0xB5]     = 0x000D9040
  1638. '    sqrt_table_even[0xB6]     = 0x000D98E8
  1639. '    sqrt_table_even[0xB7]     = 0x000DA18D
  1640. '    sqrt_table_even[0xB8]     = 0x000DAA2F
  1641. '    sqrt_table_even[0xB9]     = 0x000DB2CF
  1642. '    sqrt_table_even[0xBA]     = 0x000DBB6D
  1643. '    sqrt_table_even[0xBB]     = 0x000DC408
  1644. '    sqrt_table_even[0xBC]     = 0x000DCCA0
  1645. '    sqrt_table_even[0xBD]     = 0x000DD536
  1646. '    sqrt_table_even[0xBE]     = 0x000DDDCA
  1647. '    sqrt_table_even[0xBF]     = 0x000DE65B
  1648. '    sqrt_table_even[0xC0]     = 0x000DEEEA
  1649. '    sqrt_table_even[0xC1]     = 0x000DF776
  1650. '    sqrt_table_even[0xC2]     = 0x000E0000
  1651. '    sqrt_table_even[0xC3]     = 0x000E0887
  1652. '    sqrt_table_even[0xC4]     = 0x000E110C
  1653. '    sqrt_table_even[0xC5]     = 0x000E198E
  1654. '    sqrt_table_even[0xC6]     = 0x000E220E
  1655. '    sqrt_table_even[0xC7]     = 0x000E2A8C
  1656. '    sqrt_table_even[0xC8]     = 0x000E3307
  1657. '    sqrt_table_even[0xC9]     = 0x000E3B80
  1658. '    sqrt_table_even[0xCA]     = 0x000E43F7
  1659. '    sqrt_table_even[0xCB]     = 0x000E4C6B
  1660. '    sqrt_table_even[0xCC]     = 0x000E54DD
  1661. '    sqrt_table_even[0xCD]     = 0x000E5D4C
  1662. '    sqrt_table_even[0xCE]     = 0x000E65B9
  1663. '    sqrt_table_even[0xCF]     = 0x000E6E24
  1664. '    sqrt_table_even[0xD0]     = 0x000E768D
  1665. '    sqrt_table_even[0xD1]     = 0x000E7EF3
  1666. '    sqrt_table_even[0xD2]     = 0x000E8757
  1667. '    sqrt_table_even[0xD3]     = 0x000E8FB8
  1668. '    sqrt_table_even[0xD4]     = 0x000E9818
  1669. '    sqrt_table_even[0xD5]     = 0x000EA075
  1670. '    sqrt_table_even[0xD6]     = 0x000EA8CF
  1671. '    sqrt_table_even[0xD7]     = 0x000EB128
  1672. '    sqrt_table_even[0xD8]     = 0x000EB97E
  1673. '    sqrt_table_even[0xD9]     = 0x000EC1D2
  1674. '    sqrt_table_even[0xDA]     = 0x000ECA23
  1675. '    sqrt_table_even[0xDB]     = 0x000ED273
  1676. '    sqrt_table_even[0xDC]     = 0x000EDAC0
  1677. '    sqrt_table_even[0xDD]     = 0x000EE30B
  1678. '    sqrt_table_even[0xDE]     = 0x000EEB53
  1679. '    sqrt_table_even[0xDF]     = 0x000EF39A
  1680. '    sqrt_table_even[0xE0]     = 0x000EFBDE
  1681. '    sqrt_table_even[0xE1]     = 0x000F0420
  1682. '    sqrt_table_even[0xE2]     = 0x000F0C60
  1683. '    sqrt_table_even[0xE3]     = 0x000F149E
  1684. '    sqrt_table_even[0xE4]     = 0x000F1CD9
  1685. '    sqrt_table_even[0xE5]     = 0x000F2513
  1686. '    sqrt_table_even[0xE6]     = 0x000F2D4A
  1687. '    sqrt_table_even[0xE7]     = 0x000F357F
  1688. '    sqrt_table_even[0xE8]     = 0x000F3DB2
  1689. '    sqrt_table_even[0xE9]     = 0x000F45E2
  1690. '    sqrt_table_even[0xEA]     = 0x000F4E11
  1691. '    sqrt_table_even[0xEB]     = 0x000F563D
  1692. '    sqrt_table_even[0xEC]     = 0x000F5E67
  1693. '    sqrt_table_even[0xED]     = 0x000F6690
  1694. '    sqrt_table_even[0xEE]     = 0x000F6EB6
  1695. '    sqrt_table_even[0xEF]     = 0x000F76DA
  1696. '    sqrt_table_even[0xF0]     = 0x000F7EFB
  1697. '    sqrt_table_even[0xF1]     = 0x000F871B
  1698. '    sqrt_table_even[0xF2]     = 0x000F8F39
  1699. '    sqrt_table_even[0xF3]     = 0x000F9754
  1700. '    sqrt_table_even[0xF4]     = 0x000F9F6E
  1701. '    sqrt_table_even[0xF5]     = 0x000FA785
  1702. '    sqrt_table_even[0xF6]     = 0x000FAF9B
  1703. '    sqrt_table_even[0xF7]     = 0x000FB7AE
  1704. '    sqrt_table_even[0xF8]     = 0x000FBFBF
  1705. '    sqrt_table_even[0xF9]     = 0x000FC7CE
  1706. '    sqrt_table_even[0xFA]     = 0x000FCFDB
  1707. '    sqrt_table_even[0xFB]     = 0x000FD7E6
  1708. '    sqrt_table_even[0xFC]     = 0x000FDFEF
  1709. '    sqrt_table_even[0xFD]     = 0x000FE7F6
  1710. '    sqrt_table_even[0xFE]     = 0x000FEFFB
  1711. '    sqrt_table_even[0xFF]     = 0x000FF7FE
  1712. 'END SUB
  1713. 'END FUNCTION
  1714. '
  1715. '
  1716. ' ********************
  1717. ' *****  TAN ()  *****  Hart #4287
  1718. ' ********************
  1719. '
  1720. 'FUNCTION DOUBLE TAN (DOUBLE a) DOUBLE
  1721. '
  1722. ' FPU routine
  1723. '
  1724. '    k# = XxxFPTAN (a, @j#)
  1725. '    RETURN (k# / j#)
  1726. '
  1727. ' hand coded routine
  1728. '
  1729. '    fold into 0 - 360
  1730. '
  1731. '    DO WHILE (a < 0#)
  1732. '        a = a + $$TWOPI
  1733. '    LOOP
  1734. '
  1735. '    DO WHILE (a >= $$TWOPI)
  1736. '        a = a - $$TWOPI
  1737. '    LOOP
  1738. '
  1739. '    fold into 0-90 with appropriate sign
  1740. '
  1741. '    SELECT CASE TRUE
  1742. '        CASE (a <= $$PIDIV2)    : theSign = +1#
  1743. '        CASE (a <= $$PI)            : theSign = -1#    : a = $$PI - a
  1744. '        CASE (a <= $$PI3DIV2)    : theSign = +1#    : a = a - $$PI
  1745. '        CASE ELSE                            : theSign = -1#    : a = $$TWOPI - a
  1746. '    END SELECT
  1747. '
  1748. '    IF (a = $$PIDIV2) THEN RETURN ($$PINF)
  1749. '
  1750. '    a = a / $$PIDIV4
  1751. '    aa = a * a
  1752. '    x = aa * .1751083054422421995601107123d-1 - .279527948722905226792457575d2
  1753. '    x = x * aa + .6171941889092866899718078509d4
  1754. '    x = x * aa - .3498924461690337314553322577d6
  1755. '    x = x * aa + .413160917052212537541564775d7
  1756. '    y = aa - .4976002053470988434868766867d3                ' was aa * 1# - 497.xxxxx
  1757. '    y = y * aa + .5497880219885277215781502314d5
  1758. '    y = y * aa - .1527149650334576959585193088d7
  1759. '    y = y * aa + .5260528179299214050832459853d7
  1760. '    RETURN (((a * x) / y) * theSign)
  1761. 'END FUNCTION
  1762. '
  1763. '
  1764. ' *********************
  1765. ' *****  TANH ()  *****     >>>  may have problems similar to XmaSinh ???  <<<
  1766. ' *********************
  1767. '
  1768. FUNCTION DOUBLE TANH (DOUBLE v) DOUBLE
  1769.     SELECT CASE TRUE
  1770.         CASE (v = 0#)            :                                        RETURN (0#)
  1771.         CASE (v >= 19#)        :                                        RETURN (1#)
  1772.         CASE (v <= -19#)    :                                        RETURN (-1#)
  1773.         CASE (v > .1#)        : vv = v+v                : RETURN ((EXP(vv) - 1#) / (EXP(vv) + 1#))
  1774.         CASE (v < -.1#)        : vv = v+v                : RETURN ((EXP(vv) - 1#) / (EXP(vv) + 1#))
  1775.         CASE ELSE                    : xv = Expmo(v+v)    : RETURN (xv / (2# + xv))
  1776.     END SELECT
  1777. END FUNCTION
  1778. '
  1779. '
  1780. ' **********************
  1781. ' *****  Asin0 ()  *****  Hart #4731  :  INTERNAL FUNCTION  :  returns RADIANS
  1782. ' **********************
  1783. '
  1784. FUNCTION DOUBLE Asin0 (DOUBLE aa) DOUBLE
  1785. '
  1786.     x = aa * .28887940520319958009# - .1532823762188072258146d2#
  1787.     x = x * aa + .15627431676469845164879d3#
  1788.     x = x * aa - .61608581811333824880753d3#
  1789.     x = x * aa + .1131354445141400374754542d4#
  1790.     x = x * aa - .975678343527571443444814d3#
  1791.     x = x * aa + .31952141659008684896086d3#
  1792.     y = aa - .282183566472129245114d2#
  1793.     y = y * aa + .22430629957413222990564d3#
  1794.     y = y * aa - .76632677210477257595839d3#
  1795.     y = y * aa + .1278878991056988003191528d4#
  1796.     y = y * aa - .1028931912959252211748113d4#
  1797.     y = y * aa + .319521416590086848208615d3#
  1798.     RETURN (x / y)
  1799. END FUNCTION
  1800. '
  1801. '
  1802. ' **********************
  1803. ' *****  Atan0 ()  *****  Hart #5034 : after scaling. called from XmaAtan
  1804. ' **********************
  1805. '
  1806. '        INTERNAL FUNCTION--returns RADIANS
  1807. '
  1808. FUNCTION DOUBLE Atan0 (DOUBLE t) DOUBLE
  1809.   tt = t * t
  1810.   x = tt * .2070905893536336532# + .48914801854996690693d1
  1811.   x = x * tt + .16006919587592563754615d2
  1812.   x = x * tt + .125696450645933755049118d2
  1813. '
  1814.   y = tt + .91098182645306962582d1                            ' was tt * .1d1 + .91xxxx
  1815.   y = y * tt + .20196801275790307967877d2
  1816.   y = y * tt + .125696450645933755237905d2
  1817.   RETURN ((t * x) / y)
  1818. END FUNCTION
  1819. '
  1820. '
  1821. ' *********************
  1822. ' *****  Exp0 ()  *****  Hart #1324
  1823. ' *********************
  1824. '
  1825. FUNCTION DOUBLE Exp0 (DOUBLE v) DOUBLE
  1826.   vv = v * v
  1827.     x = vv * .606133079074800425748489607d2 + .3028561978211645920624269927d5
  1828.     x = x * vv + .2080283036505962712855955242d7
  1829.     y = vv + .1749220769510571455899141717d4
  1830.     y = y * vv + .3277095471932811805340200719d6
  1831.     y = y * vv + .6002428040825173665336946908d7
  1832.     z = (((2# * (v * x)) / (y - (v * x))) + 1#)
  1833.     RETURN (z)
  1834. '
  1835. '
  1836. '    Hart #1067        (original Exp0)
  1837. '    vv = v * v
  1838. '    x = vv * .23093347753750233624d-1 + .20202065651286927227886d2
  1839. '    x = x * vv + .1513906799054338915894328d4
  1840. '    y = vv + .233184211427481623790295d3
  1841. '    y = y * vv + .4368211662727558498496814d4
  1842. '    z = (y + (v * x)) / (y - (v * x))
  1843. '
  1844. '
  1845. '    Hart #1069        (original Exp1)
  1846. '    vv = v * v
  1847. '    x = vv * .6061485330061080841615584556d2 + .3028697169744036299076048876d5
  1848. '    x = x * vv + .2080384346694663001443843411d7
  1849. '    y = vv + .1749287689093076403844945335d4
  1850. '    y = y * vv + .3277251518082914423057964422d6
  1851. '    y = y * vv + .6002720360238832528230907598d7
  1852. '    z = (y + (v * x)) / (y - (v * x))
  1853. END FUNCTION
  1854. '
  1855. '
  1856. ' **********************
  1857. ' *****  Expmo ()  *****  Hart #1802
  1858. ' **********************
  1859. '
  1860. FUNCTION DOUBLE Expmo (DOUBLE v) DOUBLE
  1861.     vv = v * v
  1862.     x = vv * .3333206802628149222225154d-1 + .1400022777377555304975874306d2
  1863.     x = x * vv + .5040127554054843027460596926d3
  1864.     y = vv + .1120025814484651564577585494d3
  1865.     y = y * vv + .1008025510810968605492139581d4
  1866.     z = (2# * v * x) / (y - (v * x))
  1867.     RETURN (z)
  1868. END FUNCTION
  1869. '
  1870. '
  1871. ' *********************
  1872. ' *****  Log0 ()  *****  Hart #2705
  1873. ' *********************
  1874. '
  1875. FUNCTION DOUBLE Log0 (DOUBLE v) DOUBLE
  1876.   vv = v * v
  1877.     x = vv * .4210873712179797145# - .96376909336868659324d1
  1878.     x = x * vv + .30957292821537650062264d2
  1879.     x = x * vv - .240139179559210509868484d2
  1880.     y = vv - .89111090279378312337d1
  1881.     y = y * vv + .19480966070088973051623d2
  1882.     y = y * vv - .120069589779605254717525d2
  1883.     z = x / y
  1884.     RETURN (z)
  1885. '
  1886. '
  1887. '    Hart #2665  (original Log0)
  1888. '
  1889. '    vv    = v * v
  1890. '    x        = vv * .16948212488d0 + .1811136267967d0
  1891. '    x        = x * vv + .22223823332791d0
  1892. '    x        = x * vv + .2857140915904889d0
  1893. '    x        = x * vv + .400000001206045365d0
  1894. '    x        = x * vv + .6666666666633660894d0
  1895. '    x        = x * vv + .200000000000000261007d1
  1896. END FUNCTION
  1897. '
  1898. '
  1899. ' *********************
  1900. ' *****  Clog ()  *****
  1901. ' *********************
  1902. '
  1903. FUNCTION DOUBLE Clog (DOUBLE v) DOUBLE
  1904.     SLONG upper, lower, exp, exp0, exp1
  1905.     AUTOX SLONG cexp
  1906.     XLONG  ##ERROR
  1907. '
  1908.     K1 = DMAKE (0xBF2BD010, 0x5C610CA9)
  1909.     K2 = SMAKE (0x3F318000)
  1910. '
  1911.     SELECT CASE TRUE
  1912.         CASE (v <= 0#)    : ##ERROR = $$ErrorNatureInvalidArgument : RETURN (0#)
  1913.         CASE (v  = 1#)    : RETURN (0#)
  1914.     END SELECT
  1915. '
  1916.     upper    = DHIGH (v)
  1917.     lower    = DLOW (v)
  1918.     exp        = upper {11, 20}
  1919.     exp0    = upper AND 0x800FFFFF
  1920.     exp1    = exp0 OR 0x3FE00000
  1921.     exp2    = exp - 1022
  1922.     frac    = DMAKE (exp1, lower)                        ' frac is between 1/2 and 1
  1923.  
  1924. '    fric  = frexp(v, &cexp)
  1925. '    PRINT "frac, exp2 = ";; HEX$(DHIGH(frac), 8);; HEX$(DLOW(frac), 8);;; exp2
  1926. '    PRINT "fric, cexp = ";; HEX$(DHIGH(fric), 8);; HEX$(DLOW(fric), 8);;; cexp
  1927. '
  1928.     IF (frac >= $$INVSQRT2) THEN                    ' frac must be between 1/sqrt2 and sqrt2
  1929.         z = (frac - 1#) / (frac + 1#)                ' it is!
  1930.     ELSE
  1931.         DEC exp2                                                        ' nope: dec the exponent, transfer to frac
  1932.         z = (frac - .5#) / (frac + .5#)            ' mult frac by 2 (clever...)
  1933.     END IF
  1934.     zp = Clog0 (z)
  1935. '    PRINT "Clog:  exp2 * ln2, zp = "; exp2 * $$LOGE2;; zp
  1936.     j = ((exp2 * K1) + zp) + (exp2 * K2)
  1937. '    j = exp2 * $$LOGE2 + zp
  1938.     RETURN (j)
  1939. END FUNCTION
  1940. '
  1941. '
  1942. ' **********************
  1943. ' *****  Clog0 ()  *****
  1944. ' **********************
  1945. '
  1946. FUNCTION DOUBLE Clog0 (DOUBLE x) DOUBLE
  1947. '
  1948.     N1 = DMAKE (0xBFE94415, 0xB356BD28)
  1949.     N2 = DMAKE (0x4030624A, 0x2016AFED)
  1950.     N3 = DMAKE (0xC05007FF, 0x12B3B59B)
  1951.     D1 = DMAKE (0x3FF00000, 0x00000000)
  1952.     D2 = DMAKE (0xC041D580, 0x4B67CE0E)
  1953.     D3 = DMAKE (0x40738083, 0xFA15267E)
  1954.     D4 = DMAKE (0xC0880BFE, 0x9C0D9077)
  1955. '
  1956.     v        = 2# * x
  1957.     vv    = v * v
  1958.     vvv = v * vv
  1959.     x        = (((vv * N1 + N2) * vv) + N3) * vvv
  1960.     y   = (((((vv * D1) + D2) * vv) + D3) * vv) + D4
  1961.     z   = (x / y) + v
  1962.     RETURN (z)
  1963. END FUNCTION
  1964. END PROGRAM
  1965.