home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2239 < prev    next >
Encoding:
Internet Message Format  |  1990-12-28  |  41.2 KB
  1. From: glenn@extro.ucc.su.OZ.AU (Glenn Geers)
  2. Newsgroups: alt.sources
  3. Subject: 80386 alternative math library part02/04
  4. Message-ID: <1990Dec7.215654.495@metro.ucc.su.OZ.AU>
  5. Date: 7 Dec 90 21:56:54 GMT
  6.  
  7. Submitted-by: glenn@trantor
  8. Archive-name: libfpu/part02
  9.  
  10. #!/bin/sh
  11. # This is part 02 of libfpu
  12. # ============= fabs.s ==============
  13. if test -f 'fabs.s' -a X"$1" != X"-c"; then
  14.     echo 'x - skipping fabs.s (File already exists)'
  15. else
  16. echo 'x - extracting fabs.s (Text)'
  17. sed 's/^X//' << 'SHAR_EOF' > 'fabs.s' &&
  18. /*
  19. ** This file is part of the alternative 80386 math library and is
  20. ** covered by the GNU General Public license with my modification
  21. ** as noted in the README file that accompanied this file.
  22. **
  23. ** Copyright 1990 G. Geers
  24. **
  25. */
  26. X
  27. X    .align 4
  28. X    .globl fabs
  29. fabs:
  30. X    pushl %ebp
  31. X    movl %esp,%ebp
  32. X
  33. X    fldl 8(%ebp)
  34. X    fabs
  35. X
  36. X    leave
  37. X    ret
  38. SHAR_EOF
  39. chmod 0644 fabs.s ||
  40. echo 'restore of fabs.s failed'
  41. Wc_c="`wc -c < 'fabs.s'`"
  42. test 322 -eq "$Wc_c" ||
  43.     echo 'fabs.s: original size 322, current size' "$Wc_c"
  44. fi
  45. # ============= finite.s ==============
  46. if test -f 'finite.s' -a X"$1" != X"-c"; then
  47.     echo 'x - skipping finite.s (File already exists)'
  48. else
  49. echo 'x - extracting finite.s (Text)'
  50. sed 's/^X//' << 'SHAR_EOF' > 'finite.s' &&
  51. /*
  52. ** This file is part of the alternative 80386 math library and is
  53. ** covered by the GNU General Public license with my modification
  54. ** as noted in the README file that accompanied this file.
  55. **
  56. ** Copyright 1990 G. Geers
  57. **
  58. */
  59. X
  60. X    .align 4
  61. X    .globl finite
  62. finite:
  63. X    pushl %ebp
  64. X    movl %esp,%ebp
  65. X
  66. X    movl 12(%ebp), %eax
  67. X    andl $0x7ff00000, %eax
  68. X    cmpl $0x7ff00000, %eax
  69. X    je    .Lnotfinite
  70. X
  71. X    movl $1, %eax
  72. X    leave
  73. X    ret
  74. X
  75. .Lnotfinite:
  76. X    movl $0, %eax
  77. X    leave
  78. X    ret
  79. SHAR_EOF
  80. chmod 0644 finite.s ||
  81. echo 'restore of finite.s failed'
  82. Wc_c="`wc -c < 'finite.s'`"
  83. test 447 -eq "$Wc_c" ||
  84.     echo 'finite.s: original size 447, current size' "$Wc_c"
  85. fi
  86. # ============= floor.s ==============
  87. if test -f 'floor.s' -a X"$1" != X"-c"; then
  88.     echo 'x - skipping floor.s (File already exists)'
  89. else
  90. echo 'x - extracting floor.s (Text)'
  91. sed 's/^X//' << 'SHAR_EOF' > 'floor.s' &&
  92. /*
  93. ** This file is part of the alternative 80386 math library and is
  94. ** covered by the GNU General Public license with my modification
  95. ** as noted in the README file that accompanied this file.
  96. **
  97. ** Copyright 1990 G. Geers
  98. **
  99. */
  100. X
  101. X    .align 4
  102. X    .globl floor
  103. floor:
  104. X    pushl %ebp
  105. X    movl %esp,%ebp
  106. X    subl $8, %esp
  107. X    
  108. X    fldl 8(%ebp)        /* load data */
  109. X
  110. X    fstcw -12(%ebp)        /* store control word */
  111. X    fstcw -16(%ebp)        /* store it again */
  112. X    orw $0x0400, -16(%ebp)    /* round toward -inf */
  113. X    fldcw -16(%ebp)        /* store new control word */
  114. X
  115. X    frndint            /* rounding gives floor(x) */
  116. X    fldcw -12(%ebp)        /* restore original control word */
  117. X
  118. X    leave
  119. X    ret
  120. SHAR_EOF
  121. chmod 0644 floor.s ||
  122. echo 'restore of floor.s failed'
  123. Wc_c="`wc -c < 'floor.s'`"
  124. test 628 -eq "$Wc_c" ||
  125.     echo 'floor.s: original size 628, current size' "$Wc_c"
  126. fi
  127. # ============= fmod.s ==============
  128. if test -f 'fmod.s' -a X"$1" != X"-c"; then
  129.     echo 'x - skipping fmod.s (File already exists)'
  130. else
  131. echo 'x - extracting fmod.s (Text)'
  132. sed 's/^X//' << 'SHAR_EOF' > 'fmod.s' &&
  133. /*
  134. ** This file is part of the alternative 80386 math library and is
  135. ** covered by the GNU General Public license with my modification
  136. ** as noted in the README file that accompanied this file.
  137. **
  138. ** Copyright 1990 G. Geers
  139. **
  140. */
  141. X
  142. X    .align 4
  143. X    .globl fmod
  144. fmod:
  145. X    pushl %ebp
  146. X    movl %esp,%ebp
  147. X    
  148. X    fldl 16(%ebp)
  149. X    fldl 8(%ebp)
  150. .Lnotred:
  151. X    fprem
  152. X
  153. X    fstsw %ax
  154. X    sahf
  155. X    jp .Lnotred
  156. X
  157. X    leave
  158. X    ret
  159. SHAR_EOF
  160. chmod 0644 fmod.s ||
  161. echo 'restore of fmod.s failed'
  162. Wc_c="`wc -c < 'fmod.s'`"
  163. test 380 -eq "$Wc_c" ||
  164.     echo 'fmod.s: original size 380, current size' "$Wc_c"
  165. fi
  166. # ============= fpumath.h ==============
  167. if test -f 'fpumath.h' -a X"$1" != X"-c"; then
  168.     echo 'x - skipping fpumath.h (File already exists)'
  169. else
  170. echo 'x - extracting fpumath.h (Text)'
  171. sed 's/^X//' << 'SHAR_EOF' > 'fpumath.h' &&
  172. /*
  173. ** This file is covered by the GNU General Public Licence
  174. ** 
  175. ** This file should be installed somewhere in your standard include
  176. ** search path and given the name fpumath.h.
  177. */
  178. X
  179. /*
  180. ** This is a good way of flagging whether 
  181. ** you are using the alternate stuff
  182. */
  183. #define __FPU__
  184. X
  185. extern double j0(), j1(), jn(), y0(), y1(), yn();
  186. extern double erf(), erfc();
  187. extern double exp(), log(), log10(), pow(), sqrt();
  188. extern double expm1(), log1p(), exp10();
  189. extern double exp2(), log2();
  190. extern double floor(), ceil(), fabs();
  191. extern double lgamma(), gamma();
  192. extern double sinh(), cosh(), tanh();
  193. extern double asinh(), acosh(), atanh();
  194. extern double sin(), cos(), tan(), asin();
  195. extern double acos(), atan(), atan2(), hypot();
  196. extern double sqrtp();
  197. X
  198. /*
  199. ** IEEE functions
  200. */
  201. X
  202. extern double rint(), drem(), scalb(), logb(), copysign();
  203. extern double infinity(), nextafter();
  204. extern int isnan(), iszero(), isinf(), isnormal(), issubnormal();
  205. extern int signbit(), finite();
  206. extern void ieee_retrospective();
  207. extern double max_normal(), min_normal(), min_subnormal(), max_subnormal();
  208. extern double quiet_nan(),signaling_nan();
  209. extern double nextafter();
  210. X
  211. /*
  212. ** Extensions
  213. */
  214. X
  215. extern int setinternal(), setcont();
  216. X
  217. X
  218. /*
  219. ** Useful defines for setcont
  220. */
  221. X
  222. #define MASK_ALL 0x127f
  223. #define MASK_ALL1 0x137f     /* extra precision */
  224. X
  225. #define __infinity infinity
  226. X
  227. #define M_E    2.7182818284590452354
  228. #define M_LOG2E    1.4426950408889634074
  229. #define M_LOG10E    0.43429448190325182765
  230. #define M_LN2    0.69314718055994530942
  231. #define M_LN10    2.30258509299404568402
  232. #define M_PI    3.14159265358979323846
  233. #define M_PI_2    1.57079632679489661923
  234. #define M_PI_4    0.78539816339744830962
  235. #define M_1_PI    0.31830988618379067154
  236. #define M_2_PI    0.63661977236758134308
  237. #define M_2_SQRTPI    1.12837916709551257390
  238. #define M_SQRT2    1.41421356237309504880
  239. #define M_SQRT1_2    0.70710678118654752440
  240. X
  241. #define MAXFLOAT    ((float)3.40282346638528860e+38)
  242. #define    MINFLOAT    ((float)1.40129846432481707e-45)
  243. X
  244. #ifndef IEEE
  245. #define    MAXDOUBLE    1.79769313486231570e+308  /* wrong in math.h */
  246. #define    MINDOUBLE    4.94065645841246544e-324
  247. #else
  248. #define    MAXDOUBLE    (max_normal())
  249. #define    MINDOUBLE    (min_subnormal())
  250. #endif
  251. X
  252. #define ULP        1.1102230246251565e-16
  253. X
  254. #define    HUGE_VAL    3.40282346638528860e+38
  255. X
  256. #ifdef IEEE
  257. #define HUGE    (__infinity())
  258. #else
  259. #define HUGE     MAXDOUBLE
  260. #endif
  261. SHAR_EOF
  262. chmod 0644 fpumath.h ||
  263. echo 'restore of fpumath.h failed'
  264. Wc_c="`wc -c < 'fpumath.h'`"
  265. test 2314 -eq "$Wc_c" ||
  266.     echo 'fpumath.h: original size 2314, current size' "$Wc_c"
  267. fi
  268. # ============= gamma.c ==============
  269. if test -f 'gamma.c' -a X"$1" != X"-c"; then
  270.     echo 'x - skipping gamma.c (File already exists)'
  271. else
  272. echo 'x - extracting gamma.c (Text)'
  273. sed 's/^X//' << 'SHAR_EOF' > 'gamma.c' &&
  274. /*
  275. ** This file is part of the alternative 80386 math library and is
  276. ** covered by the GNU General Public license with my modification
  277. ** as noted in the README file that accompanied this file.
  278. **
  279. ** This implementation takes no notice of the fact that gamma(z) is undefined
  280. ** if z is 0 or a negative integer - beware.
  281. **
  282. ** Copyright 1990 G. Geers
  283. **
  284. */
  285. X
  286. /*
  287. ** The limits are approximate
  288. */
  289. #define U_LIM 100.0
  290. #define L_LIM -100.0
  291. X
  292. extern double lgamma();
  293. extern double exp();
  294. extern int signgam;
  295. X
  296. double gamma(x)
  297. double x;
  298. {
  299. X    if (x >= L_LIM && x <= U_LIM)
  300. X        return((double)signgam*exp(lgamma(x)));
  301. }
  302. SHAR_EOF
  303. chmod 0644 gamma.c ||
  304. echo 'restore of gamma.c failed'
  305. Wc_c="`wc -c < 'gamma.c'`"
  306. test 604 -eq "$Wc_c" ||
  307.     echo 'gamma.c: original size 604, current size' "$Wc_c"
  308. fi
  309. # ============= hypot.s ==============
  310. if test -f 'hypot.s' -a X"$1" != X"-c"; then
  311.     echo 'x - skipping hypot.s (File already exists)'
  312. else
  313. echo 'x - extracting hypot.s (Text)'
  314. sed 's/^X//' << 'SHAR_EOF' > 'hypot.s' &&
  315. /*
  316. ** This file is part of the alternative 80386 math library and is
  317. ** covered by the GNU General Public license with my modification
  318. ** as noted in the README file that accompanied this file.
  319. **
  320. ** Copyright 1990 G. Geers
  321. **
  322. */
  323. X
  324. X    .align 4
  325. .globl hypot
  326. hypot:
  327. X    pushl %ebp
  328. X    movl %esp,%ebp
  329. X
  330. X    fldl 8(%ebp)
  331. X    fmull 8(%ebp)
  332. X    fldl 16(%ebp)
  333. X    fmull 16(%ebp)
  334. X    faddp
  335. X    fsqrt
  336. X
  337. X    leave
  338. X    ret
  339. SHAR_EOF
  340. chmod 0644 hypot.s ||
  341. echo 'restore of hypot.s failed'
  342. Wc_c="`wc -c < 'hypot.s'`"
  343. test 377 -eq "$Wc_c" ||
  344.     echo 'hypot.s: original size 377, current size' "$Wc_c"
  345. fi
  346. # ============= ieee_ext.s ==============
  347. if test -f 'ieee_ext.s' -a X"$1" != X"-c"; then
  348.     echo 'x - skipping ieee_ext.s (File already exists)'
  349. else
  350. echo 'x - extracting ieee_ext.s (Text)'
  351. sed 's/^X//' << 'SHAR_EOF' > 'ieee_ext.s' &&
  352. /*
  353. ** This file is part of the alternative 80386 math library and is
  354. ** covered by the GNU General Public license with my modification
  355. ** as noted in the README file that accompanied this file.
  356. **
  357. ** Copyright 1990 G. Geers
  358. **
  359. */
  360. X
  361. X    .align 4
  362. X    .globl isnan
  363. isnan:
  364. X    pushl %ebp
  365. X    movl %esp,%ebp
  366. X
  367. X    movl 12(%ebp), %eax
  368. X    andl $0x7ff00000, %eax
  369. X    cmpl $0x7ff00000, %eax
  370. X    jne .Lnotnan
  371. X    movl 12(%ebp), %eax
  372. X    andl $0xfffff, %eax
  373. X    orl 8(%ebp), %eax
  374. X    je .Lnotnan
  375. X
  376. X    movl $1, %eax
  377. X    leave 
  378. X    ret
  379. X
  380. .Lnotnan:
  381. X    movl $0, %eax
  382. X
  383. .Ldone:
  384. X    leave
  385. X    ret
  386. X
  387. X    .align 4
  388. X    .globl isinf
  389. isinf:
  390. X    pushl %ebp
  391. X    movl %esp,%ebp
  392. X
  393. X    movl 12(%ebp), %eax
  394. X    andl $0x7ff00000, %eax
  395. X    cmpl $0x7ff00000, %eax
  396. X    je .Lcouldbeinf
  397. X
  398. .Lnotinf:
  399. X    movl $0, %eax
  400. X    leave
  401. X    ret
  402. X
  403. .Lcouldbeinf:
  404. X    andl $0xfffff, %eax
  405. X    orl 8(%ebp), %eax
  406. X    jne .Lnotinf
  407. X
  408. X    movl $1, %eax
  409. X    leave
  410. X    ret
  411. X
  412. X    .align 4
  413. X    .globl iszero
  414. iszero:
  415. X    pushl %ebp
  416. X    movl %esp,%ebp
  417. X
  418. X    movl 12(%ebp), %eax
  419. X    cmpl $0x0, %eax
  420. X    je .Lcouldbezero
  421. .Lnotzero:
  422. X    movl $0, %eax
  423. X    leave
  424. X    ret
  425. X
  426. .Lcouldbezero:
  427. X    andl $0xfffff, %eax
  428. X    orl 8(%ebp), %eax
  429. X    jne .Lnotzero
  430. X
  431. X    movl $1, %eax
  432. X    leave
  433. X    ret
  434. X
  435. X    .align 4
  436. X    .globl signbit
  437. signbit:
  438. X    pushl %ebp
  439. X    movl %esp,%ebp
  440. X
  441. X    movl 12(%ebp), %eax
  442. X    andl $0x80000000, %eax
  443. X    cmpl $0x80000000, %eax
  444. X    jne .Lpos
  445. X    movl $1, %eax
  446. X    leave
  447. X    ret
  448. X
  449. .Lpos:
  450. X    movl $0, %eax
  451. X    leave
  452. X    ret
  453. X
  454. X    .align 4
  455. X    .globl issubnormal
  456. issubnormal:
  457. X    pushl %ebp
  458. X    movl %esp,%ebp
  459. X
  460. X    movl 12(%ebp), %eax
  461. X    andl $0x7ff00000, %eax
  462. X    cmpl $0x0, %eax
  463. X    je .Lcouldbesub
  464. X
  465. .Lnotsubnorm:
  466. X    movl $0, %eax
  467. X    leave
  468. X    ret
  469. X
  470. .Lcouldbesub:
  471. X    movl 12(%ebp), %eax
  472. X    andl $0xfffff, %eax
  473. X    orl 8(%ebp), %eax
  474. X    je .Lnotsubnorm
  475. X
  476. X    movl $1, %eax
  477. X    leave
  478. X    ret
  479. X
  480. X    .align 4
  481. X    .globl isnormal
  482. isnormal:
  483. X    pushl %ebp
  484. X    movl %esp,%ebp
  485. X
  486. X    movl 12(%ebp), %eax
  487. X    andl $0x7ff00000, %eax /* mask sign bit */
  488. X    xorl $0x7ff00000, %eax 
  489. X    cmpl $0x0, %eax
  490. X    je .Lnotnorm
  491. X    cmpl $0x7ff00000, %eax
  492. X    je .Lcouldbenorm
  493. X
  494. .Lnorm:
  495. X    movl $1, %eax
  496. X    leave
  497. X    ret
  498. X
  499. .Lcouldbenorm:
  500. X    movl 12(%ebp), %eax
  501. X    andl $0xfffff, %eax
  502. X    orl 8(%ebp), %eax
  503. X    je .Lnorm
  504. X
  505. .Lnotnorm:
  506. X    movl $0, %eax
  507. X    leave
  508. X    ret
  509. SHAR_EOF
  510. chmod 0644 ieee_ext.s ||
  511. echo 'restore of ieee_ext.s failed'
  512. Wc_c="`wc -c < 'ieee_ext.s'`"
  513. test 1977 -eq "$Wc_c" ||
  514.     echo 'ieee_ext.s: original size 1977, current size' "$Wc_c"
  515. fi
  516. # ============= ieee_retro.c ==============
  517. if test -f 'ieee_retro.c' -a X"$1" != X"-c"; then
  518.     echo 'x - skipping ieee_retro.c (File already exists)'
  519. else
  520. echo 'x - extracting ieee_retro.c (Text)'
  521. sed 's/^X//' << 'SHAR_EOF' > 'ieee_retro.c' &&
  522. /*
  523. ** This file is part of the alternative 80386 math library and is
  524. ** covered by the GNU General Public license with my modification
  525. ** as noted in the README file that accompanied this file.
  526. **
  527. ** Copyright 1990 G. Geers
  528. **
  529. */
  530. X
  531. #include <stdio.h>
  532. X
  533. ieee_retrospective(f)
  534. FILE *f;
  535. {
  536. X    int status;
  537. X
  538. X    if (f == (FILE *)NULL)
  539. X        f = stderr;
  540. X
  541. X    status = _getsw();
  542. X
  543. X    if (status & 0xdf) {
  544. X        fprintf(f,"\n");
  545. X        fprintf(f,"Warning: the following IEEE floating-point");
  546. X        fprintf(f," arithmetic exceptions\n");
  547. X        fprintf(f,"occurred in this program and were never cleared:\n");
  548. X        if (status & 0x0001) fprintf(f,"Invalid Operation;\n");
  549. X        if (status & 0x0002) fprintf(f,"Denormal;\n");
  550. X        if (status & 0x0004) fprintf(f,"Division by Zero;\n");
  551. X        if (status & 0x0008) fprintf(f,"Overflow;\n");
  552. X        if (status & 0x0010) fprintf(f,"Underflow;\n");
  553. X        fprintf(f,"\n");
  554. X    }
  555. }
  556. SHAR_EOF
  557. chmod 0644 ieee_retro.c ||
  558. echo 'restore of ieee_retro.c failed'
  559. Wc_c="`wc -c < 'ieee_retro.c'`"
  560. test 853 -eq "$Wc_c" ||
  561.     echo 'ieee_retro.c: original size 853, current size' "$Wc_c"
  562. fi
  563. # ============= ieee_values.s ==============
  564. if test -f 'ieee_values.s' -a X"$1" != X"-c"; then
  565.     echo 'x - skipping ieee_values.s (File already exists)'
  566. else
  567. echo 'x - extracting ieee_values.s (Text)'
  568. sed 's/^X//' << 'SHAR_EOF' > 'ieee_values.s' &&
  569. /*
  570. ** This file is part of the alternative 80386 math library and is
  571. ** covered by the GNU General Public license with my modification
  572. ** as noted in the README file that accompanied this file.
  573. **
  574. ** Copyright 1990 G. Geers
  575. **
  576. */
  577. X
  578. X    .align 4
  579. X    .globl max_normal
  580. X
  581. max_normal:
  582. X    pushl %ebp
  583. X    movl %esp,%ebp
  584. X    subl $8, %esp
  585. X
  586. X    movl $0x7fefffff, -12(%ebp)
  587. X    movl $0xffffffff, -16(%ebp)
  588. X
  589. X    fldl -16(%ebp)
  590. X    
  591. X    leave
  592. X    ret
  593. X
  594. X    .align 4
  595. X    .globl min_normal
  596. X
  597. min_normal:
  598. X    pushl %ebp
  599. X    movl %esp,%ebp
  600. X    subl $8, %esp
  601. X
  602. X    movl $0x00100000, -12(%ebp)
  603. X    movl $0x00000001, -16(%ebp)
  604. X
  605. X    fldl -16(%ebp)
  606. X    
  607. X    leave
  608. X    ret
  609. X
  610. X    .align 4
  611. X    .globl min_subnormal
  612. X
  613. min_subnormal:
  614. X    pushl %ebp
  615. X    movl %esp,%ebp
  616. X    subl $8, %esp
  617. X
  618. X    movl $0x0, -12(%ebp)
  619. X    movl $0x00000001, -16(%ebp)
  620. X
  621. X    fldl -16(%ebp)
  622. X    
  623. X    leave
  624. X    ret
  625. X
  626. X    .align 4
  627. X    .globl max_subnormal
  628. X
  629. max_subnormal:
  630. X    pushl %ebp
  631. X    movl %esp,%ebp
  632. X    subl $8, %esp
  633. X
  634. X    movl $0x000fffff, -12(%ebp)
  635. X    movl $0xffffffff, -16(%ebp)
  636. X
  637. X    fldl -16(%ebp)
  638. X    
  639. X    leave
  640. X    ret
  641. X
  642. X    .align 4
  643. X    .globl quiet_nan
  644. X
  645. quiet_nan:
  646. X    pushl %ebp
  647. X    movl %esp,%ebp
  648. X    subl $8, %esp
  649. X
  650. X    movl $0x7fffffff, -12(%ebp)
  651. X    movl $0xffffffff, -16(%ebp)
  652. X
  653. X    fldl -16(%ebp)
  654. X    
  655. X    leave
  656. X    ret
  657. X
  658. X    .align 4
  659. X    .globl signaling_nan
  660. X
  661. signaling_nan:
  662. X    pushl %ebp
  663. X    movl %esp,%ebp
  664. X    subl $8, %esp
  665. X
  666. X    movl $0x7ff00000, -12(%ebp)
  667. X    movl $0x00000001, -16(%ebp)
  668. X
  669. X    fldl -16(%ebp)
  670. X    
  671. X    leave
  672. X    ret
  673. SHAR_EOF
  674. chmod 0644 ieee_values.s ||
  675. echo 'restore of ieee_values.s failed'
  676. Wc_c="`wc -c < 'ieee_values.s'`"
  677. test 1295 -eq "$Wc_c" ||
  678.     echo 'ieee_values.s: original size 1295, current size' "$Wc_c"
  679. fi
  680. # ============= infinity.s ==============
  681. if test -f 'infinity.s' -a X"$1" != X"-c"; then
  682.     echo 'x - skipping infinity.s (File already exists)'
  683. else
  684. echo 'x - extracting infinity.s (Text)'
  685. sed 's/^X//' << 'SHAR_EOF' > 'infinity.s' &&
  686. /*
  687. ** This file is part of the alternative 80386 math library and is
  688. ** covered by the GNU General Public license with my modification
  689. ** as noted in the README file that accompanied this file.
  690. **
  691. ** Copyright 1990 G. Geers
  692. **
  693. */
  694. X
  695. X    .align 4
  696. X    .globl infinity
  697. infinity:
  698. X    pushl %ebp
  699. X    movl %esp,%ebp
  700. X    subl $8, %esp
  701. X
  702. X    movl $0x7ff00000, -12(%ebp)
  703. X    movl $0x0, -16(%ebp)
  704. X
  705. X    fldl -16(%ebp)
  706. X    
  707. X    leave
  708. X    ret
  709. SHAR_EOF
  710. chmod 0644 infinity.s ||
  711. echo 'restore of infinity.s failed'
  712. Wc_c="`wc -c < 'infinity.s'`"
  713. test 394 -eq "$Wc_c" ||
  714.     echo 'infinity.s: original size 394, current size' "$Wc_c"
  715. fi
  716. # ============= j0.c ==============
  717. if test -f 'j0.c' -a X"$1" != X"-c"; then
  718.     echo 'x - skipping j0.c (File already exists)'
  719. else
  720. echo 'x - extracting j0.c (Text)'
  721. sed 's/^X//' << 'SHAR_EOF' > 'j0.c' &&
  722. /*
  723. X * Copyright (c) 1985 Regents of the University of California.
  724. X * All rights reserved.  The Berkeley software License Agreement
  725. X * specifies the terms and conditions for redistribution.
  726. X */
  727. X
  728. #ifndef lint
  729. static char sccsid[] = "@(#)j0.c    5.3 (Berkeley) 9/22/88";
  730. #endif /* not lint */
  731. X
  732. /*
  733. X    floating point Bessel's function
  734. X    of the first and second kinds
  735. X    of order zero
  736. X
  737. X    j0(x) returns the value of J0(x)
  738. X    for all real values of x.
  739. X
  740. X    There are no error returns.
  741. X    Calls sin, cos, sqrt.
  742. X
  743. X    There is a niggling bug in J0 which
  744. X    causes errors up to 2e-16 for x in the
  745. X    interval [-8,8].
  746. X    The bug is caused by an inappropriate order
  747. X    of summation of the series.  rhm will fix it
  748. X    someday.
  749. X
  750. X    Coefficients are from Hart & Cheney.
  751. X    #5849 (19.22D)
  752. X    #6549 (19.25D)
  753. X    #6949 (19.41D)
  754. X
  755. X    y0(x) returns the value of Y0(x)
  756. X    for positive real values of x.
  757. X    For x<=0, if on the VAX, error number EDOM is set and
  758. X    the reserved operand fault is generated;
  759. X    otherwise (an IEEE machine) an invalid operation is performed.
  760. X
  761. X    Calls sin, cos, sqrt, log, j0.
  762. X
  763. X    The values of Y0 have not been checked
  764. X    to more than ten places.
  765. X
  766. X    Coefficients are from Hart & Cheney.
  767. X    #6245 (18.78D)
  768. X    #6549 (19.25D)
  769. X    #6949 (19.41D)
  770. */
  771. X
  772. #include "mathimpl.h"
  773. X
  774. extern double sin();
  775. extern double cos();
  776. extern double sqrt();
  777. extern double log();
  778. X
  779. double j0();
  780. X
  781. #if defined(vax)||defined(tahoe)
  782. #include <errno.h>
  783. #else    /* defined(vax)||defined(tahoe) */
  784. static const double zero = 0.e0;
  785. #endif    /* defined(vax)||defined(tahoe) */
  786. X
  787. static double pzero, qzero;
  788. X
  789. static const double tpi    = .6366197723675813430755350535e0;
  790. static const double pio4    = .7853981633974483096156608458e0;
  791. static const double p1[] = {
  792. X    0.4933787251794133561816813446e21,
  793. X    -.1179157629107610536038440800e21,
  794. X    0.6382059341072356562289432465e19,
  795. X    -.1367620353088171386865416609e18,
  796. X    0.1434354939140344111664316553e16,
  797. X    -.8085222034853793871199468171e13,
  798. X    0.2507158285536881945555156435e11,
  799. X    -.4050412371833132706360663322e8,
  800. X    0.2685786856980014981415848441e5,
  801. };
  802. static const double q1[] = {
  803. X    0.4933787251794133562113278438e21,
  804. X    0.5428918384092285160200195092e19,
  805. X    0.3024635616709462698627330784e17,
  806. X    0.1127756739679798507056031594e15,
  807. X    0.3123043114941213172572469442e12,
  808. X    0.6699987672982239671814028660e9,
  809. X    0.1114636098462985378182402543e7,
  810. X    0.1363063652328970604442810507e4,
  811. X    1.0
  812. };
  813. static const double p2[] = {
  814. X    0.5393485083869438325262122897e7,
  815. X    0.1233238476817638145232406055e8,
  816. X    0.8413041456550439208464315611e7,
  817. X    0.2016135283049983642487182349e7,
  818. X    0.1539826532623911470917825993e6,
  819. X    0.2485271928957404011288128951e4,
  820. X    0.0,
  821. };
  822. static const double q2[] = {
  823. X    0.5393485083869438325560444960e7,
  824. X    0.1233831022786324960844856182e8,
  825. X    0.8426449050629797331554404810e7,
  826. X    0.2025066801570134013891035236e7,
  827. X    0.1560017276940030940592769933e6,
  828. X    0.2615700736920839685159081813e4,
  829. X    1.0,
  830. };
  831. static const double p3[] = {
  832. X    -.3984617357595222463506790588e4,
  833. X    -.1038141698748464093880530341e5,
  834. X    -.8239066313485606568803548860e4,
  835. X    -.2365956170779108192723612816e4,
  836. X    -.2262630641933704113967255053e3,
  837. X    -.4887199395841261531199129300e1,
  838. X    0.0,
  839. };
  840. static const double q3[] = {
  841. X    0.2550155108860942382983170882e6,
  842. X    0.6667454239319826986004038103e6,
  843. X    0.5332913634216897168722255057e6,
  844. X    0.1560213206679291652539287109e6,
  845. X    0.1570489191515395519392882766e5,
  846. X    0.4087714673983499223402830260e3,
  847. X    1.0,
  848. };
  849. static const double p4[] = {
  850. X    -.2750286678629109583701933175e20,
  851. X    0.6587473275719554925999402049e20,
  852. X    -.5247065581112764941297350814e19,
  853. X    0.1375624316399344078571335453e18,
  854. X    -.1648605817185729473122082537e16,
  855. X    0.1025520859686394284509167421e14,
  856. X    -.3436371222979040378171030138e11,
  857. X    0.5915213465686889654273830069e8,
  858. X    -.4137035497933148554125235152e5,
  859. };
  860. static const double q4[] = {
  861. X    0.3726458838986165881989980e21,
  862. X    0.4192417043410839973904769661e19,
  863. X    0.2392883043499781857439356652e17,
  864. X    0.9162038034075185262489147968e14,
  865. X    0.2613065755041081249568482092e12,
  866. X    0.5795122640700729537480087915e9,
  867. X    0.1001702641288906265666651753e7,
  868. X    0.1282452772478993804176329391e4,
  869. X    1.0,
  870. };
  871. X
  872. static void asympt();
  873. X
  874. double
  875. j0(arg) double arg;{
  876. X    double argsq, n, d;
  877. X    int i;
  878. X
  879. X    if(arg < 0.) arg = -arg;
  880. X    if(arg > 8.){
  881. X        asympt(arg);
  882. X        n = arg - pio4;
  883. X        return(sqrt(tpi/arg)*(pzero*cos(n) - qzero*sin(n)));
  884. X    }
  885. X    argsq = arg*arg;
  886. X    for(n=0,d=0,i=8;i>=0;i--){
  887. X        n = n*argsq + p1[i];
  888. X        d = d*argsq + q1[i];
  889. X    }
  890. X    return(n/d);
  891. }
  892. X
  893. double
  894. y0(arg) double arg;{
  895. X    double argsq, n, d;
  896. X    int i;
  897. X
  898. X    if(arg <= 0.){
  899. #if defined(vax)||defined(tahoe)
  900. X        return(infnan(EDOM));        /* NaN */
  901. #else    /* defined(vax)||defined(tahoe) */
  902. X        return(zero/zero);    /* IEEE machines: invalid operation */
  903. #endif    /* defined(vax)||defined(tahoe) */
  904. X    }
  905. X    if(arg > 8.){
  906. X        asympt(arg);
  907. X        n = arg - pio4;
  908. X        return(sqrt(tpi/arg)*(pzero*sin(n) + qzero*cos(n)));
  909. X    }
  910. X    argsq = arg*arg;
  911. X    for(n=0,d=0,i=8;i>=0;i--){
  912. X        n = n*argsq + p4[i];
  913. X        d = d*argsq + q4[i];
  914. X    }
  915. X    return(n/d + tpi*j0(arg)*log(arg));
  916. }
  917. X
  918. static void
  919. asympt(arg) double arg;{
  920. X    double zsq, n, d;
  921. X    int i;
  922. X    zsq = 64./(arg*arg);
  923. X    for(n=0,d=0,i=6;i>=0;i--){
  924. X        n = n*zsq + p2[i];
  925. X        d = d*zsq + q2[i];
  926. X    }
  927. X    pzero = n/d;
  928. X    for(n=0,d=0,i=6;i>=0;i--){
  929. X        n = n*zsq + p3[i];
  930. X        d = d*zsq + q3[i];
  931. X    }
  932. X    qzero = (8./arg)*(n/d);
  933. }
  934. SHAR_EOF
  935. chmod 0644 j0.c ||
  936. echo 'restore of j0.c failed'
  937. Wc_c="`wc -c < 'j0.c'`"
  938. test 5099 -eq "$Wc_c" ||
  939.     echo 'j0.c: original size 5099, current size' "$Wc_c"
  940. fi
  941. # ============= j1.c ==============
  942. if test -f 'j1.c' -a X"$1" != X"-c"; then
  943.     echo 'x - skipping j1.c (File already exists)'
  944. else
  945. echo 'x - extracting j1.c (Text)'
  946. sed 's/^X//' << 'SHAR_EOF' > 'j1.c' &&
  947. /*
  948. X * Copyright (c) 1985 Regents of the University of California.
  949. X * All rights reserved.  The Berkeley software License Agreement
  950. X * specifies the terms and conditions for redistribution.
  951. X */
  952. X
  953. #ifndef lint
  954. static char sccsid[] = "@(#)j1.c    5.3 (Berkeley) 9/22/88";
  955. #endif /* not lint */
  956. X
  957. /*
  958. X    floating point Bessel's function
  959. X    of the first and second kinds
  960. X    of order one
  961. X
  962. X    j1(x) returns the value of J1(x)
  963. X    for all real values of x.
  964. X
  965. X    There are no error returns.
  966. X    Calls sin, cos, sqrt.
  967. X
  968. X    There is a niggling bug in J1 which
  969. X    causes errors up to 2e-16 for x in the
  970. X    interval [-8,8].
  971. X    The bug is caused by an inappropriate order
  972. X    of summation of the series.  rhm will fix it
  973. X    someday.
  974. X
  975. X    Coefficients are from Hart & Cheney.
  976. X    #6050 (20.98D)
  977. X    #6750 (19.19D)
  978. X    #7150 (19.35D)
  979. X
  980. X    y1(x) returns the value of Y1(x)
  981. X    for positive real values of x.
  982. X    For x<=0, if on the VAX, error number EDOM is set and
  983. X    the reserved operand fault is generated;
  984. X    otherwise (an IEEE machine) an invalid operation is performed.
  985. X
  986. X    Calls sin, cos, sqrt, log, j1.
  987. X
  988. X    The values of Y1 have not been checked
  989. X    to more than ten places.
  990. X
  991. X    Coefficients are from Hart & Cheney.
  992. X    #6447 (22.18D)
  993. X    #6750 (19.19D)
  994. X    #7150 (19.35D)
  995. */
  996. X
  997. #include "mathimpl.h"
  998. X
  999. extern double sin();
  1000. extern double cos();
  1001. extern double sqrt();
  1002. extern double log();
  1003. X
  1004. double j1();
  1005. X
  1006. #if defined(vax)||defined(tahoe)
  1007. #include <errno.h>
  1008. #else    /* defined(vax)||defined(tahoe) */
  1009. static const double zero = 0.e0;
  1010. #endif    /* defined(vax)||defined(tahoe) */
  1011. X
  1012. static double pzero, qzero;
  1013. X
  1014. static const double tpi    = .6366197723675813430755350535e0;
  1015. static const double pio4    = .7853981633974483096156608458e0;
  1016. static const double p1[] = {
  1017. X    0.581199354001606143928050809e21,
  1018. X    -.6672106568924916298020941484e20,
  1019. X    0.2316433580634002297931815435e19,
  1020. X    -.3588817569910106050743641413e17,
  1021. X    0.2908795263834775409737601689e15,
  1022. X    -.1322983480332126453125473247e13,
  1023. X    0.3413234182301700539091292655e10,
  1024. X    -.4695753530642995859767162166e7,
  1025. X    0.2701122710892323414856790990e4,
  1026. };
  1027. static const double q1[] = {
  1028. X    0.1162398708003212287858529400e22,
  1029. X    0.1185770712190320999837113348e20,
  1030. X    0.6092061398917521746105196863e17,
  1031. X    0.2081661221307607351240184229e15,
  1032. X    0.5243710262167649715406728642e12,
  1033. X    0.1013863514358673989967045588e10,
  1034. X    0.1501793594998585505921097578e7,
  1035. X    0.1606931573481487801970916749e4,
  1036. X    1.0,
  1037. };
  1038. static const double p2[] = {
  1039. X    -.4435757816794127857114720794e7,
  1040. X    -.9942246505077641195658377899e7,
  1041. X    -.6603373248364939109255245434e7,
  1042. X    -.1523529351181137383255105722e7,
  1043. X    -.1098240554345934672737413139e6,
  1044. X    -.1611616644324610116477412898e4,
  1045. X    0.0,
  1046. };
  1047. static const double q2[] = {
  1048. X    -.4435757816794127856828016962e7,
  1049. X    -.9934124389934585658967556309e7,
  1050. X    -.6585339479723087072826915069e7,
  1051. X    -.1511809506634160881644546358e7,
  1052. X    -.1072638599110382011903063867e6,
  1053. X    -.1455009440190496182453565068e4,
  1054. X    1.0,
  1055. };
  1056. static const double p3[] = {
  1057. X    0.3322091340985722351859704442e5,
  1058. X    0.8514516067533570196555001171e5,
  1059. X    0.6617883658127083517939992166e5,
  1060. X    0.1849426287322386679652009819e5,
  1061. X    0.1706375429020768002061283546e4,
  1062. X    0.3526513384663603218592175580e2,
  1063. X    0.0,
  1064. };
  1065. static const double q3[] = {
  1066. X    0.7087128194102874357377502472e6,
  1067. X    0.1819458042243997298924553839e7,
  1068. X    0.1419460669603720892855755253e7,
  1069. X    0.4002944358226697511708610813e6,
  1070. X    0.3789022974577220264142952256e5,
  1071. X    0.8638367769604990967475517183e3,
  1072. X    1.0,
  1073. };
  1074. static const double p4[] = {
  1075. X    -.9963753424306922225996744354e23,
  1076. X    0.2655473831434854326894248968e23,
  1077. X    -.1212297555414509577913561535e22,
  1078. X    0.2193107339917797592111427556e20,
  1079. X    -.1965887462722140658820322248e18,
  1080. X    0.9569930239921683481121552788e15,
  1081. X    -.2580681702194450950541426399e13,
  1082. X    0.3639488548124002058278999428e10,
  1083. X    -.2108847540133123652824139923e7,
  1084. X    0.0,
  1085. };
  1086. static const double q4[] = {
  1087. X    0.5082067366941243245314424152e24,
  1088. X    0.5435310377188854170800653097e22,
  1089. X    0.2954987935897148674290758119e20,
  1090. X    0.1082258259408819552553850180e18,
  1091. X    0.2976632125647276729292742282e15,
  1092. X    0.6465340881265275571961681500e12,
  1093. X    0.1128686837169442121732366891e10,
  1094. X    0.1563282754899580604737366452e7,
  1095. X    0.1612361029677000859332072312e4,
  1096. X    1.0,
  1097. };
  1098. X
  1099. static void asympt();
  1100. X
  1101. double
  1102. j1(arg) double arg;{
  1103. X    double xsq, n, d, x;
  1104. X    int i;
  1105. X
  1106. X    x = arg;
  1107. X    if(x < 0.) x = -x;
  1108. X    if(x > 8.){
  1109. X        asympt(x);
  1110. X        n = x - 3.*pio4;
  1111. X        n = sqrt(tpi/x)*(pzero*cos(n) - qzero*sin(n));
  1112. X        if(arg <0.) n = -n;
  1113. X        return(n);
  1114. X    }
  1115. X    xsq = x*x;
  1116. X    for(n=0,d=0,i=8;i>=0;i--){
  1117. X        n = n*xsq + p1[i];
  1118. X        d = d*xsq + q1[i];
  1119. X    }
  1120. X    return(arg*n/d);
  1121. }
  1122. X
  1123. double
  1124. y1(arg) double arg;{
  1125. X    double xsq, n, d, x;
  1126. X    int i;
  1127. X
  1128. X    x = arg;
  1129. X    if(x <= 0.){
  1130. #if defined(vax)||defined(tahoe)
  1131. X        return(infnan(EDOM));        /* NaN */
  1132. #else    /* defined(vax)||defined(tahoe) */
  1133. X        return(zero/zero);    /* IEEE machines: invalid operation */
  1134. #endif    /* defined(vax)||defined(tahoe) */
  1135. X    }
  1136. X    if(x > 8.){
  1137. X        asympt(x);
  1138. X        n = x - 3*pio4;
  1139. X        return(sqrt(tpi/x)*(pzero*sin(n) + qzero*cos(n)));
  1140. X    }
  1141. X    xsq = x*x;
  1142. X    for(n=0,d=0,i=9;i>=0;i--){
  1143. X        n = n*xsq + p4[i];
  1144. X        d = d*xsq + q4[i];
  1145. X    }
  1146. X    return(x*n/d + tpi*(j1(x)*log(x)-1./x));
  1147. }
  1148. X
  1149. static void
  1150. asympt(arg) double arg;{
  1151. X    double zsq, n, d;
  1152. X    int i;
  1153. X    zsq = 64./(arg*arg);
  1154. X    for(n=0,d=0,i=6;i>=0;i--){
  1155. X        n = n*zsq + p2[i];
  1156. X        d = d*zsq + q2[i];
  1157. X    }
  1158. X    pzero = n/d;
  1159. X    for(n=0,d=0,i=6;i>=0;i--){
  1160. X        n = n*zsq + p3[i];
  1161. X        d = d*zsq + q3[i];
  1162. X    }
  1163. X    qzero = (8./arg)*(n/d);
  1164. }
  1165. SHAR_EOF
  1166. chmod 0644 j1.c ||
  1167. echo 'restore of j1.c failed'
  1168. Wc_c="`wc -c < 'j1.c'`"
  1169. test 5169 -eq "$Wc_c" ||
  1170.     echo 'j1.c: original size 5169, current size' "$Wc_c"
  1171. fi
  1172. # ============= jn.c ==============
  1173. if test -f 'jn.c' -a X"$1" != X"-c"; then
  1174.     echo 'x - skipping jn.c (File already exists)'
  1175. else
  1176. echo 'x - extracting jn.c (Text)'
  1177. sed 's/^X//' << 'SHAR_EOF' > 'jn.c' &&
  1178. /*
  1179. X * Copyright (c) 1985 Regents of the University of California.
  1180. X * All rights reserved.  The Berkeley software License Agreement
  1181. X * specifies the terms and conditions for redistribution.
  1182. X */
  1183. X
  1184. #ifndef lint
  1185. static char sccsid[] = "@(#)jn.c    5.3 (Berkeley) 9/22/88";
  1186. #endif /* not lint */
  1187. X
  1188. /*
  1189. X    floating point Bessel's function of
  1190. X    the first and second kinds and of
  1191. X    integer order.
  1192. X
  1193. X    int n;
  1194. X    double x;
  1195. X    jn(n,x);
  1196. X
  1197. X    returns the value of Jn(x) for all
  1198. X    integer values of n and all real values
  1199. X    of x.
  1200. X
  1201. X    There are no error returns.
  1202. X    Calls j0, j1.
  1203. X
  1204. X    For n=0, j0(x) is called,
  1205. X    for n=1, j1(x) is called,
  1206. X    for n<x, forward recursion us used starting
  1207. X    from values of j0(x) and j1(x).
  1208. X    for n>x, a continued fraction approximation to
  1209. X    j(n,x)/j(n-1,x) is evaluated and then backward
  1210. X    recursion is used starting from a supposed value
  1211. X    for j(n,x). The resulting value of j(0,x) is
  1212. X    compared with the actual value to correct the
  1213. X    supposed value of j(n,x).
  1214. X
  1215. X    yn(n,x) is similar in all respects, except
  1216. X    that forward recursion is used for all
  1217. X    values of n>1.
  1218. */
  1219. X
  1220. #include "mathimpl.h"
  1221. X
  1222. extern double j0();
  1223. extern double j1();
  1224. X
  1225. #if defined(vax)||defined(tahoe)
  1226. #include <errno.h>
  1227. #else    /* defined(vax)||defined(tahoe) */
  1228. static double zero = 0.e0;
  1229. #endif    /* defined(vax)||defined(tahoe) */
  1230. X
  1231. double
  1232. jn(n,x) int n; double x;{
  1233. X    int i;
  1234. X    double a, b, temp;
  1235. X    double xsq, t;
  1236. X
  1237. X    if(n<0){
  1238. X        n = -n;
  1239. X        x = -x;
  1240. X    }
  1241. X    if(n==0) return(j0(x));
  1242. X    if(n==1) return(j1(x));
  1243. X    if(x == 0.) return(0.);
  1244. X    if(n>x) goto recurs;
  1245. X
  1246. X    a = j0(x);
  1247. X    b = j1(x);
  1248. X    for(i=1;i<n;i++){
  1249. X        temp = b;
  1250. X        b = (2.*i/x)*b - a;
  1251. X        a = temp;
  1252. X    }
  1253. X    return(b);
  1254. X
  1255. recurs:
  1256. X    xsq = x*x;
  1257. X    for(t=0,i=n+16;i>n;i--){
  1258. X        t = xsq/(2.*i - t);
  1259. X    }
  1260. X    t = x/(2.*n-t);
  1261. X
  1262. X    a = t;
  1263. X    b = 1;
  1264. X    for(i=n-1;i>0;i--){
  1265. X        temp = b;
  1266. X        b = (2.*i/x)*b - a;
  1267. X        a = temp;
  1268. X    }
  1269. X    return(t*j0(x)/b);
  1270. }
  1271. X
  1272. double
  1273. yn(n,x) int n; double x;{
  1274. X    int i;
  1275. X    int sign;
  1276. X    double a, b, temp;
  1277. X
  1278. X    if (x <= 0) {
  1279. #if defined(vax)||defined(tahoe)
  1280. X        return(infnan(EDOM));    /* NaN */
  1281. #else    /* defined(vax)||defined(tahoe) */
  1282. X        return(zero/zero);    /* IEEE machines: invalid operation */
  1283. #endif    /* defined(vax)||defined(tahoe) */
  1284. X    }
  1285. X    sign = 1;
  1286. X    if(n<0){
  1287. X        n = -n;
  1288. X        if(n%2 == 1) sign = -1;
  1289. X    }
  1290. X    if(n==0) return(y0(x));
  1291. X    if(n==1) return(sign*y1(x));
  1292. X
  1293. X    a = y0(x);
  1294. X    b = y1(x);
  1295. X    for(i=1;i<n;i++){
  1296. X        temp = b;
  1297. X        b = (2.*i/x)*b - a;
  1298. X        a = temp;
  1299. X    }
  1300. X    return(sign*b);
  1301. }
  1302. SHAR_EOF
  1303. chmod 0644 jn.c ||
  1304. echo 'restore of jn.c failed'
  1305. Wc_c="`wc -c < 'jn.c'`"
  1306. test 2310 -eq "$Wc_c" ||
  1307.     echo 'jn.c: original size 2310, current size' "$Wc_c"
  1308. fi
  1309. # ============= lgamma.c ==============
  1310. if test -f 'lgamma.c' -a X"$1" != X"-c"; then
  1311.     echo 'x - skipping lgamma.c (File already exists)'
  1312. else
  1313. echo 'x - extracting lgamma.c (Text)'
  1314. sed 's/^X//' << 'SHAR_EOF' > 'lgamma.c' &&
  1315. /*
  1316. X * Copyright (c) 1985 Regents of the University of California.
  1317. X * All rights reserved.  The Berkeley software License Agreement
  1318. X * specifies the terms and conditions for redistribution.
  1319. X */
  1320. X
  1321. #ifndef lint
  1322. static char sccsid[] = "@(#)lgamma.c    5.3 (Berkeley) 9/22/88";
  1323. #endif /* not lint */
  1324. X
  1325. /*
  1326. X    C program for floating point log Gamma function
  1327. X
  1328. X    lgamma(x) computes the log of the absolute
  1329. X    value of the Gamma function.
  1330. X    The sign of the Gamma function is returned in the
  1331. X    external quantity signgam.
  1332. X
  1333. X    The coefficients for expansion around zero
  1334. X    are #5243 from Hart & Cheney; for expansion
  1335. X    around infinity they are #5404.
  1336. X
  1337. X    Calls log, floor and sin.
  1338. */
  1339. X
  1340. #include "mathimpl.h"
  1341. X
  1342. extern double log();
  1343. extern double floor();
  1344. extern double sin();
  1345. X
  1346. #if defined(vax)||defined(tahoe)
  1347. #include <errno.h>
  1348. #endif    /* defined(vax)||defined(tahoe) */
  1349. X
  1350. int    signgam = 0;
  1351. static const double goobie = 0.9189385332046727417803297;  /* log(2*pi)/2 */
  1352. static const double pi       = 3.1415926535897932384626434;
  1353. X
  1354. #define M 6
  1355. #define N 8
  1356. static const double p1[] = {
  1357. X    0.83333333333333101837e-1,
  1358. X    -.277777777735865004e-2,
  1359. X    0.793650576493454e-3,
  1360. X    -.5951896861197e-3,
  1361. X    0.83645878922e-3,
  1362. X    -.1633436431e-2,
  1363. };
  1364. static const double p2[] = {
  1365. X    -.42353689509744089647e5,
  1366. X    -.20886861789269887364e5,
  1367. X    -.87627102978521489560e4,
  1368. X    -.20085274013072791214e4,
  1369. X    -.43933044406002567613e3,
  1370. X    -.50108693752970953015e2,
  1371. X    -.67449507245925289918e1,
  1372. X    0.0,
  1373. };
  1374. static const double q2[] = {
  1375. X    -.42353689509744090010e5,
  1376. X    -.29803853309256649932e4,
  1377. X    0.99403074150827709015e4,
  1378. X    -.15286072737795220248e4,
  1379. X    -.49902852662143904834e3,
  1380. X    0.18949823415702801641e3,
  1381. X    -.23081551524580124562e2,
  1382. X    0.10000000000000000000e1,
  1383. };
  1384. X
  1385. static double pos(), neg(), asym();
  1386. X
  1387. double
  1388. lgamma(arg)
  1389. double arg;
  1390. {
  1391. X
  1392. X    signgam = 1.;
  1393. X    if(arg <= 0.) return(neg(arg));
  1394. X    if(arg > 8.) return(asym(arg));
  1395. X    return(log(pos(arg)));
  1396. }
  1397. X
  1398. static double
  1399. asym(arg)
  1400. double arg;
  1401. {
  1402. X    double n, argsq;
  1403. X    int i;
  1404. X
  1405. X    argsq = 1./(arg*arg);
  1406. X    for(n=0,i=M-1; i>=0; i--){
  1407. X        n = n*argsq + p1[i];
  1408. X    }
  1409. X    return((arg-.5)*log(arg) - arg + goobie + n/arg);
  1410. }
  1411. X
  1412. static double
  1413. neg(arg)
  1414. double arg;
  1415. {
  1416. X    double t;
  1417. X
  1418. X    arg = -arg;
  1419. X     /*
  1420. X      * to see if arg were a true integer, the old code used the
  1421. X      * mathematically correct observation:
  1422. X      * sin(n*pi) = 0 <=> n is an integer.
  1423. X      * but in finite precision arithmetic, sin(n*PI) will NEVER
  1424. X      * be zero simply because n*PI is a rational number.  hence
  1425. X      *    it failed to work with our newer, more accurate sin()
  1426. X      * which uses true pi to do the argument reduction...
  1427. X      *    temp = sin(pi*arg);
  1428. X      */
  1429. X    t = floor(arg);
  1430. X    if (arg - t  > 0.5e0)
  1431. X        t += 1.e0;                /* t := integer nearest arg */
  1432. #if defined(vax)||defined(tahoe)
  1433. X    if (arg == t) {
  1434. X        return(infnan(ERANGE));        /* +INF */
  1435. X    }
  1436. #endif    /* defined(vax)||defined(tahoe) */
  1437. X    signgam = (int) (t - 2*floor(t/2));    /* signgam =  1 if t was odd, */
  1438. X                        /*            0 if t was even */
  1439. X    signgam = signgam - 1 + signgam;    /* signgam =  1 if t was odd, */
  1440. X                        /*           -1 if t was even */
  1441. X    t = arg - t;                /*  -0.5 <= t <= 0.5 */
  1442. X    if (t < 0.e0) {
  1443. X        t = -t;
  1444. X        signgam = -signgam;
  1445. X    }
  1446. X    return(-log(arg*pos(arg)*sin(pi*t)/pi));
  1447. }
  1448. X
  1449. static double
  1450. pos(arg)
  1451. double arg;
  1452. {
  1453. X    double n, d, s;
  1454. X    register i;
  1455. X
  1456. X    if(arg < 2.) return(pos(arg+1.)/arg);
  1457. X    if(arg > 3.) return((arg-1.)*pos(arg-1.));
  1458. X
  1459. X    s = arg - 2.;
  1460. X    for(n=0,d=0,i=N-1; i>=0; i--){
  1461. X        n = n*s + p2[i];
  1462. X        d = d*s + q2[i];
  1463. X    }
  1464. X    return(n/d);
  1465. }
  1466. SHAR_EOF
  1467. chmod 0644 lgamma.c ||
  1468. echo 'restore of lgamma.c failed'
  1469. Wc_c="`wc -c < 'lgamma.c'`"
  1470. test 3382 -eq "$Wc_c" ||
  1471.     echo 'lgamma.c: original size 3382, current size' "$Wc_c"
  1472. fi
  1473. # ============= log.s ==============
  1474. if test -f 'log.s' -a X"$1" != X"-c"; then
  1475.     echo 'x - skipping log.s (File already exists)'
  1476. else
  1477. echo 'x - extracting log.s (Text)'
  1478. sed 's/^X//' << 'SHAR_EOF' > 'log.s' &&
  1479. /*
  1480. ** This file is part of the alternative 80386 math library and is
  1481. ** covered by the GNU General Public license with my modification
  1482. ** as noted in the README file that accompanied this file.
  1483. **
  1484. ** Copyright 1990 G. Geers
  1485. **
  1486. */
  1487. X
  1488. X    .align 4
  1489. X    .globl log
  1490. log:
  1491. X    pushl %ebp
  1492. X    movl %esp,%ebp
  1493. X
  1494. X    fldln2
  1495. X    fldl 8(%ebp)
  1496. X    fyl2x
  1497. X
  1498. X    leave
  1499. X    ret
  1500. SHAR_EOF
  1501. chmod 0644 log.s ||
  1502. echo 'restore of log.s failed'
  1503. Wc_c="`wc -c < 'log.s'`"
  1504. test 329 -eq "$Wc_c" ||
  1505.     echo 'log.s: original size 329, current size' "$Wc_c"
  1506. fi
  1507. # ============= log10.s ==============
  1508. if test -f 'log10.s' -a X"$1" != X"-c"; then
  1509.     echo 'x - skipping log10.s (File already exists)'
  1510. else
  1511. echo 'x - extracting log10.s (Text)'
  1512. sed 's/^X//' << 'SHAR_EOF' > 'log10.s' &&
  1513. /*
  1514. ** This file is part of the alternative 80386 math library and is
  1515. ** covered by the GNU General Public license with my modification
  1516. ** as noted in the README file that accompanied this file.
  1517. **
  1518. ** Copyright 1990 G. Geers
  1519. **
  1520. */
  1521. X
  1522. X    .align 4
  1523. X    .globl log10
  1524. log10:
  1525. X    pushl %ebp
  1526. X    movl %esp,%ebp
  1527. X
  1528. X    fldlg2
  1529. X    fldl 8(%ebp)
  1530. X    fyl2x
  1531. X
  1532. X    leave
  1533. X    ret
  1534. SHAR_EOF
  1535. chmod 0644 log10.s ||
  1536. echo 'restore of log10.s failed'
  1537. Wc_c="`wc -c < 'log10.s'`"
  1538. test 333 -eq "$Wc_c" ||
  1539.     echo 'log10.s: original size 333, current size' "$Wc_c"
  1540. fi
  1541. # ============= log1p.s ==============
  1542. if test -f 'log1p.s' -a X"$1" != X"-c"; then
  1543.     echo 'x - skipping log1p.s (File already exists)'
  1544. else
  1545. echo 'x - extracting log1p.s (Text)'
  1546. sed 's/^X//' << 'SHAR_EOF' > 'log1p.s' &&
  1547. /*
  1548. ** This file is part of the alternative 80386 math library and is
  1549. ** covered by the GNU General Public license with my modification
  1550. ** as noted in the README file that accompanied this file.
  1551. **
  1552. ** Copyright 1990 G. Geers
  1553. **
  1554. */
  1555. X
  1556. X    .align 4
  1557. X    .globl log1p
  1558. log1p:
  1559. X    pushl %ebp
  1560. X    movl %esp,%ebp
  1561. X
  1562. X    fldln2
  1563. X    fldl 8(%ebp)
  1564. X    fld1
  1565. X    faddp
  1566. X    fyl2x
  1567. X
  1568. X    leave
  1569. X    ret
  1570. SHAR_EOF
  1571. chmod 0644 log1p.s ||
  1572. echo 'restore of log1p.s failed'
  1573. Wc_c="`wc -c < 'log1p.s'`"
  1574. test 346 -eq "$Wc_c" ||
  1575.     echo 'log1p.s: original size 346, current size' "$Wc_c"
  1576. fi
  1577. # ============= log2.s ==============
  1578. if test -f 'log2.s' -a X"$1" != X"-c"; then
  1579.     echo 'x - skipping log2.s (File already exists)'
  1580. else
  1581. echo 'x - extracting log2.s (Text)'
  1582. sed 's/^X//' << 'SHAR_EOF' > 'log2.s' &&
  1583. /*
  1584. ** This file is part of the alternative 80386 math library and is
  1585. ** covered by the GNU General Public license with my modification
  1586. ** as noted in the README file that accompanied this file.
  1587. **
  1588. ** Copyright 1990 G. Geers
  1589. **
  1590. */
  1591. X
  1592. X    .align 4
  1593. X    .globl log2
  1594. log2:
  1595. X    pushl %ebp
  1596. X    movl %esp,%ebp
  1597. X
  1598. X    fld1
  1599. X    fldl 8(%ebp)
  1600. X    fyl2x
  1601. X
  1602. X    leave
  1603. X    ret
  1604. SHAR_EOF
  1605. chmod 0644 log2.s ||
  1606. echo 'restore of log2.s failed'
  1607. Wc_c="`wc -c < 'log2.s'`"
  1608. test 329 -eq "$Wc_c" ||
  1609.     echo 'log2.s: original size 329, current size' "$Wc_c"
  1610. fi
  1611. # ============= logb.s ==============
  1612. if test -f 'logb.s' -a X"$1" != X"-c"; then
  1613.     echo 'x - skipping logb.s (File already exists)'
  1614. else
  1615. echo 'x - extracting logb.s (Text)'
  1616. sed 's/^X//' << 'SHAR_EOF' > 'logb.s' &&
  1617. /*
  1618. ** This file is part of the alternative 80386 math library and is
  1619. ** covered by the GNU General Public license with my modification
  1620. ** as noted in the README file that accompanied this file.
  1621. **
  1622. ** Copyright 1990 G. Geers
  1623. **
  1624. */
  1625. X
  1626. X    .align 4
  1627. X    .globl logb
  1628. logb:
  1629. X    pushl %ebp
  1630. X    movl %esp,%ebp
  1631. X    
  1632. X    fldl 8(%ebp)
  1633. X    fxtract
  1634. X    fldl %st(1)
  1635. X
  1636. X    leave
  1637. X    ret
  1638. SHAR_EOF
  1639. chmod 0644 logb.s ||
  1640. echo 'restore of logb.s failed'
  1641. Wc_c="`wc -c < 'logb.s'`"
  1642. test 339 -eq "$Wc_c" ||
  1643.     echo 'logb.s: original size 339, current size' "$Wc_c"
  1644. fi
  1645. # ============= mathimpl.h ==============
  1646. if test -f 'mathimpl.h' -a X"$1" != X"-c"; then
  1647.     echo 'x - skipping mathimpl.h (File already exists)'
  1648. else
  1649. echo 'x - extracting mathimpl.h (Text)'
  1650. sed 's/^X//' << 'SHAR_EOF' > 'mathimpl.h' &&
  1651. /*
  1652. X * Copyright (c) 1988 The Regents of the University of California.
  1653. X * All rights reserved.
  1654. X *
  1655. X * Redistribution and use in source and binary forms are permitted
  1656. X * provided that: (1) source distributions retain this entire copyright
  1657. X * notice and comment, and (2) distributions including binaries display
  1658. X * the following acknowledgement:  ``This product includes software
  1659. X * developed by the University of California, Berkeley and its contributors''
  1660. X * in the documentation or other materials provided with the distribution
  1661. X * and in all advertising materials mentioning features or use of this
  1662. X * software. Neither the name of the University nor the names of its
  1663. X * contributors may be used to endorse or promote products derived
  1664. X * from this software without specific prior written permission.
  1665. X * THIS SOFTWARE IS PROVIDED ``AS IS'' AND WITHOUT ANY EXPRESS OR
  1666. X * IMPLIED WARRANTIES, INCLUDING, WITHOUT LIMITATION, THE IMPLIED
  1667. X * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
  1668. X *
  1669. X * All recipients should regard themselves as participants in an ongoing
  1670. X * research project and hence should feel obligated to report their
  1671. X * experiences (good or bad) with these elementary function codes, using
  1672. X * the sendbug(8) program, to the authors.
  1673. X *
  1674. X *    @(#)mathimpl.h    5.2 (Berkeley) 6/1/90
  1675. X */
  1676. X
  1677. #include <math.h>
  1678. X
  1679. #if defined(__STDC__) || defined(__GNUC__)
  1680. #define const const
  1681. #else
  1682. #define const /**/
  1683. #endif
  1684. X
  1685. #if defined(vax)||defined(tahoe)
  1686. X
  1687. /* Deal with different ways to concatenate in cpp */
  1688. #  if defined(__STDC__) || defined(__GNUC__)
  1689. #    define    cat3(a,b,c) a ## b ## c
  1690. #  else
  1691. #    define    cat3(a,b,c) a/**/b/**/c
  1692. #  endif
  1693. X
  1694. /* Deal with vax/tahoe byte order issues */
  1695. #  ifdef vax
  1696. #    define    cat3t(a,b,c) cat3(a,b,c)
  1697. #  else
  1698. #    define    cat3t(a,b,c) cat3(a,c,b)
  1699. #  endif
  1700. X
  1701. #  define vccast(name) (*(const double *)(cat3(name,,x)))
  1702. X
  1703. X   /*
  1704. X    * Define a constant to high precision on a Vax or Tahoe.
  1705. X    *
  1706. X    * Args are the name to define, the decimal floating point value,
  1707. X    * four 16-bit chunks of the float value in hex
  1708. X    * (because the vax and tahoe differ in float format!), the power
  1709. X    * of 2 of the hex-float exponent, and the hex-float mantissa.
  1710. X    * Most of these arguments are not used at compile time; they are
  1711. X    * used in a post-check to make sure the constants were compiled
  1712. X    * correctly.
  1713. X    *
  1714. X    * People who want to use the constant will have to do their own
  1715. X    *     #define foo vccast(foo)
  1716. X    * since CPP cannot do this for them from inside another macro (sigh).
  1717. X    * We define "vccast" if this needs doing.
  1718. X    */
  1719. #  define vc(name, value, x1,x2,x3,x4, bexp, xval) \
  1720. X    const static long cat3(name,,x)[] = {cat3t(0x,x1,x2), cat3t(0x,x3,x4)};
  1721. X
  1722. #  define ic(name, value, bexp, xval) ;
  1723. X
  1724. #else    /* vax or tahoe */
  1725. X
  1726. X   /* Hooray, we have an IEEE machine */
  1727. #  undef vccast
  1728. #  define vc(name, value, x1,x2,x3,x4, bexp, xval) ;
  1729. X
  1730. #  define ic(name, value, bexp, xval) \
  1731. X    const static double name = value;
  1732. X
  1733. #endif    /* defined(vax)||defined(tahoe) */
  1734. SHAR_EOF
  1735. chmod 0644 mathimpl.h ||
  1736. echo 'restore of mathimpl.h failed'
  1737. Wc_c="`wc -c < 'mathimpl.h'`"
  1738. test 2996 -eq "$Wc_c" ||
  1739.     echo 'mathimpl.h: original size 2996, current size' "$Wc_c"
  1740. fi
  1741. # ============= nextafter.c ==============
  1742. if test -f 'nextafter.c' -a X"$1" != X"-c"; then
  1743.     echo 'x - skipping nextafter.c (File already exists)'
  1744. else
  1745. echo 'x - extracting nextafter.c (Text)'
  1746. sed 's/^X//' << 'SHAR_EOF' > 'nextafter.c' &&
  1747. /*
  1748. ** This file is part of the alternative 80386 math library and is
  1749. ** covered by the GNU General Public license with my modification
  1750. ** as noted in the README file that accompanied this file.
  1751. **
  1752. ** Copyright 1990 G. Geers
  1753. **
  1754. ** A mix of C and assembler - well I've got the functions so I might 
  1755. ** as well use them!
  1756. **
  1757. */
  1758. X
  1759. #include "fpumath.h"
  1760. X
  1761. asm(".align 4");
  1762. asm(".Lulp:");
  1763. asm(".double 1.1102230246251565e-16");
  1764. X
  1765. asm(".align 4");
  1766. asm(".Lulpup:");
  1767. asm(".double 2.2204460492503131e-16");
  1768. X
  1769. double
  1770. nextafter(x, y)
  1771. double x, y;
  1772. {
  1773. X    asm("sub $8, %esp");
  1774. X
  1775. X    if (isnan(x) || isnan(y))
  1776. X        return(quiet_nan());
  1777. X
  1778. X    if (isinf(x))
  1779. X        if (y > x)
  1780. X            return(-max_normal());
  1781. X        else
  1782. X        if (y < x)
  1783. X            return(max_normal());
  1784. X
  1785. X    if (x == 0.0)
  1786. X        if (y > 0.0)
  1787. X            return(min_subnormal());
  1788. X        else
  1789. X            return(-min_subnormal());
  1790. X
  1791. X    if (isnormal(x)) {
  1792. X        asm("movl 12(%ebp), %eax");
  1793. X        asm("andl $0x7ff00000, %eax");
  1794. X        asm("movl %eax, -12(%ebp)");
  1795. X        asm("movl $0x0, -16(%ebp)");
  1796. X        asm("fldl -16(%ebp)");
  1797. X
  1798. X        if (fabs(x) <= 1.0) 
  1799. X            asm("fldl .Lulp");
  1800. X        else
  1801. X            asm("fldl .Lulpup");
  1802. X
  1803. X        asm("fmulp");
  1804. X
  1805. X        if (y > x) {
  1806. X            asm("fldl 8(%ebp)");
  1807. X            asm("faddp");
  1808. X            asm("leave");
  1809. X            asm("ret");
  1810. X        }
  1811. X        if (y < x) {
  1812. X            asm("fldl 8(%ebp)");
  1813. X            asm("fsubp");
  1814. X            asm("leave");
  1815. X            asm("ret");
  1816. X        }
  1817. X    }
  1818. X    else
  1819. X    if (issubnormal(x)) {
  1820. X        if (y > x) 
  1821. X            return(x + min_subnormal());
  1822. X
  1823. X        if (y < x)
  1824. X            return(x - min_subnormal());
  1825. X    }
  1826. X    
  1827. X    return(x);
  1828. }
  1829. SHAR_EOF
  1830. chmod 0644 nextafter.c ||
  1831. echo 'restore of nextafter.c failed'
  1832. Wc_c="`wc -c < 'nextafter.c'`"
  1833. test 1392 -eq "$Wc_c" ||
  1834.     echo 'nextafter.c: original size 1392, current size' "$Wc_c"
  1835. fi
  1836. true || echo 'restore of paranoia.c failed'
  1837. echo End of part 2, continue with part 3
  1838. exit 0
  1839. --
  1840. Glenn Geers                       | "So when it's over, we're back to people.
  1841. Department of Theoretical Physics |  Just to prove that human touch can have
  1842. The University of Sydney          |  no equal."
  1843. Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'
  1844.