home *** CD-ROM | disk | FTP | other *** search
/ Source Code 1992 March / Source_Code_CD-ROM_Walnut_Creek_March_1992.iso / usenet / altsrcs / 2 / 2335 < prev    next >
Internet Message Format  |  1990-12-28  |  32KB
  1. From: glenn@qed.physics.su.OZ.AU (Glenn Geers)
  2. Newsgroups: alt.sources
  3. Subject: Alternative 80386 math lib v2.0 part06/06
  4. Message-ID: <1990Dec16.210916.28774@metro.ucc.su.OZ.AU>
  5. Date: 16 Dec 90 21:09:16 GMT
  6.  
  7. Submitted-by: root@trantor
  8. Archive-name: mathlib2.0/part06
  9.  
  10. ---- Cut Here and feed the following to sh ----
  11. #!/bin/sh
  12. # this is mathlib.06 (part 6 of mathlib2.0)
  13. # do not concatenate these parts, unpack them in order with /bin/sh
  14. # file COPYING continued
  15. #
  16. if test ! -r _shar_seq_.tmp; then
  17.     echo 'Please unpack part 1 first!'
  18.     exit 1
  19. fi
  20. (read Scheck
  21.  if test "$Scheck" != 6; then
  22.     echo Please unpack part "$Scheck" next!
  23.     exit 1
  24.  else
  25.     exit 0
  26.  fi
  27. ) < _shar_seq_.tmp || exit 1
  28. if test ! -f _shar_wnt_.tmp; then
  29.     echo 'x - still skipping COPYING'
  30. else
  31. echo 'x - continuing file COPYING'
  32. sed 's/^X//' << 'SHAR_EOF' >> 'COPYING' &&
  33. X
  34. X    This program is free software; you can redistribute it and/or modify
  35. X    it under the terms of the GNU General Public License as published by
  36. X    the Free Software Foundation; either version 1, or (at your option)
  37. X    any later version.
  38. X
  39. X    This program is distributed in the hope that it will be useful,
  40. X    but WITHOUT ANY WARRANTY; without even the implied warranty of
  41. X    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  42. X    GNU General Public License for more details.
  43. X
  44. X    You should have received a copy of the GNU General Public License
  45. X    along with this program; if not, write to the Free Software
  46. X    Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  47. X
  48. Also add information on how to contact you by electronic and paper mail.
  49. X
  50. If the program is interactive, make it output a short notice like this
  51. when it starts in an interactive mode:
  52. X
  53. X    Gnomovision version 69, Copyright (C) 19xx name of author
  54. X    Gnomovision comes with ABSOLUTELY NO WARRANTY; for details type `show w'.
  55. X    This is free software, and you are welcome to redistribute it
  56. X    under certain conditions; type `show c' for details.
  57. X
  58. The hypothetical commands `show w' and `show c' should show the
  59. appropriate parts of the General Public License.  Of course, the
  60. commands you use may be called something other than `show w' and `show
  61. c'; they could even be mouse-clicks or menu items--whatever suits your
  62. program.
  63. X
  64. You should also get your employer (if you work as a programmer) or your
  65. school, if any, to sign a "copyright disclaimer" for the program, if
  66. necessary.  Here a sample; alter the names:
  67. X
  68. X  Yoyodyne, Inc., hereby disclaims all copyright interest in the
  69. X  program `Gnomovision' (a program to direct compilers to make passes
  70. X  at assemblers) written by James Hacker.
  71. X
  72. X  <signature of Ty Coon>, 1 April 1989
  73. X  Ty Coon, President of Vice
  74. X
  75. That's all there is to it!
  76. SHAR_EOF
  77. echo 'File COPYING is complete' &&
  78. chmod 0644 COPYING ||
  79. echo 'restore of COPYING failed'
  80. Wc_c="`wc -c < 'COPYING'`"
  81. test 12488 -eq "$Wc_c" ||
  82.     echo 'COPYING: original size 12488, current size' "$Wc_c"
  83. rm -f _shar_wnt_.tmp
  84. fi
  85. # ============= COPYRIGHT ==============
  86. if test -f 'COPYRIGHT' -a X"$1" != X"-c"; then
  87.     echo 'x - skipping COPYRIGHT (File already exists)'
  88.     rm -f _shar_wnt_.tmp
  89. else
  90. > _shar_wnt_.tmp
  91. echo 'x - extracting COPYRIGHT (Text)'
  92. sed 's/^X//' << 'SHAR_EOF' > 'COPYRIGHT' &&
  93. All the assembler source files and the C source files ieee_retro.c, nextafter.c
  94. and nextafterf.c are copyright 1990 to G. Geers.
  95. X
  96. The following C source files were written by W.J. Cody of Argonne National
  97. Laboratory: erf.c, erff.c, gamma.c, gammaf.c, j0.c, j0f.c, j1.c, j1f.c,
  98. lgamma.c and lgammaf.c. The float versions (names with f appended) were 
  99. derived from the originals by G. Geers. These files are used with permission.
  100. X
  101. The original BASIC version of paranoia is copyright Professor W. M. Kahan.
  102. SHAR_EOF
  103. chmod 0644 COPYRIGHT ||
  104. echo 'restore of COPYRIGHT failed'
  105. Wc_c="`wc -c < 'COPYRIGHT'`"
  106. test 504 -eq "$Wc_c" ||
  107.     echo 'COPYRIGHT: original size 504, current size' "$Wc_c"
  108. rm -f _shar_wnt_.tmp
  109. fi
  110. # ============= Makefile ==============
  111. if test -f 'Makefile' -a X"$1" != X"-c"; then
  112.     echo 'x - skipping Makefile (File already exists)'
  113.     rm -f _shar_wnt_.tmp
  114. else
  115. > _shar_wnt_.tmp
  116. echo 'x - extracting Makefile (Text)'
  117. sed 's/^X//' << 'SHAR_EOF' > 'Makefile' &&
  118. # This file is part of the 80386 alternative math library and is covered by
  119. # the GNU General Public Licence
  120. X
  121. GCC=gcc
  122. CC=cc
  123. CFLAGS=-O -fstrength-reduce -DIEEE
  124. LIBDIR=/lib
  125. INCDIR=/usr/include
  126. # if you want to use setinternal in paranoia; define PTEST
  127. #PTEST=
  128. PTEST=-DTEST
  129. VERS=2.0
  130. X
  131. .s.o:
  132. X    $(GCC) -c $<
  133. X
  134. .c.o:
  135. X    $(GCC) $(CFLAGS) -c $<
  136. X
  137. CSRC=j0.c lgamma.c gamma.c j1.c erf.c nextafter.c
  138. X
  139. CFSRC=j0f.c lgammaf.c gammaf.c j1f.c erff.c nextafterf.c
  140. X
  141. SSRC=acos.s copysign.s drem.s fabs.s hypot.s logb.s scalb.s tan.s \
  142. X    asin.s ceil.s exp.s expm1.s finite.s log.s log1p.s \
  143. X    pow.s sin.s atan.s cos.s exp10.s floor.s log10.s rint.s sqrt.s \
  144. X    exp2.s log2.s sinh.s cosh.s tanh.s \
  145. X    asinh.s acosh.s atanh.s atan2.s fmod.s ieee_ext.s \
  146. X    infinity.s sqrtp.s ieee_values.s
  147. X
  148. SFSRC=acosf.s copysignf.s fabsf.s log10f.s sinf.s acoshf.s cosf.s \
  149. X    finitef.s log1pf.s sinhf.s asinf.s coshf.s floorf.s log2f.s \
  150. X    sqrtf.s asinhf.s dremf.s fmodf.s logbf.s sqrtpf.s atan2f.s \
  151. X    exp10f.s hypotf.s logf.s tanf.s atanf.s exp2f.s ieee_extf.s \
  152. X    powf.s tanhf.s atanhf.s expf.s ieee_valuesf.s rintf.s ceilf.s \
  153. X    expm1f.s infinityf.s scalbf.s
  154. X
  155. COMMONSRC=ieee_retro.c _getsw.s setcont.s setinternal.s
  156. X
  157. ADD=CHANGELOG COPYING COPYRIGHT Makefile PROBLEMS README TODO d2dcomb.summ
  158. X
  159. ALLSRC=$(CSRC) $(CFSRC) $(SSRC) $(SFSRC) $(COMMONSRC) $(TESTSRC) $(ADD)
  160. X
  161. LIBOBJ=acos.o atanh.o erf.o hypot.o log10.o sin.o \
  162. X       acosh.o ceil.o exp.o j0.o logb.o sinh.o \
  163. X       asin.o copysign.o exp10.o expm1.o j1.o sqrt.o \
  164. X       asinh.o cos.o fabs.o pow.o tan.o \
  165. X       atan.o cosh.o finite.o lgamma.o rint.o tanh.o \
  166. X       atan2.o drem.o floor.o log.o scalb.o log1p.o gamma.o \
  167. X       exp2.o log2.o fmod.o ieee_ext.o \
  168. X       infinity.o sqrtp.o ieee_values.o nextafter.o 
  169. X
  170. LIBFOBJ=acosf.o atanhf.o erff.o hypotf.o log10f.o sinf.o \
  171. X       acoshf.o ceilf.o expf.o j0f.o logbf.o sinhf.o \
  172. X       asinf.o copysignf.o exp10f.o expm1f.o j1f.o sqrtf.o \
  173. X       asinhf.o cosf.o fabsf.o powf.o tanf.o \
  174. X       atanf.o coshf.o finitef.o lgammaf.o rintf.o tanhf.o \
  175. X       atan2f.o dremf.o floorf.o logf.o scalbf.o log1pf.o gammaf.o \
  176. X       exp2f.o log2f.o fmodf.o ieee_extf.o \
  177. X       infinityf.o sqrtpf.o ieee_valuesf.o nextafterf.o 
  178. X
  179. COMMONOBJ=ieee_retro.o _getsw.o setcont.o setinternal.o
  180. X
  181. TESTSRC=paranoia.c
  182. X
  183. all: libfpu libffpu libcfpu test ftest
  184. X
  185. libfpu: $(LIBOBJ) $(COMMONOBJ)
  186. X    ar vr libfpu.a $(LIBOBJ) $(COMMONOBJ)
  187. X
  188. libffpu: $(LIBFOBJ) $(COMMONOBJ)
  189. X    ar vr libffpu.a $(LIBFOBJ) $(COMMONOBJ)
  190. X
  191. libcfpu: $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ)
  192. X    ar vr libcfpu.a $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ)
  193. X
  194. test: libfpu $(TESTOBJ)
  195. X    $(GCC) -o paranoia $(TESTSRC) -DTEST libfpu.a
  196. X
  197. ftest: libffpu $(TESTOBJ)
  198. X    $(GCC) -o fparanoia -DTEST -DSingle $(TESTSRC) libffpu.a
  199. X
  200. paranoia: test ftest
  201. X
  202. install: libfpu libffpu libcfpu
  203. X    -cp libfpu.a $(LIBDIR)/libfpu.a
  204. X    -chown bin $(LIBDIR)/libfpu.a
  205. X    -chgrp bin $(LIBDIR)/libfpu.a
  206. X    -chmod 444 $(LIBDIR)/libfpu.a
  207. X    -cp libffpu.a $(LIBDIR)/libffpu.a
  208. X    -chown bin $(LIBDIR)/libffpu.a
  209. X    -chgrp bin $(LIBDIR)/libffpu.a
  210. X    -chmod 444 $(LIBDIR)/libffpu.a
  211. X    -cp libcfpu.a $(LIBDIR)/libcfpu.a
  212. X    -chown bin $(LIBDIR)/libcfpu.a
  213. X    -chgrp bin $(LIBDIR)/libcfpu.a
  214. X    -chmod 444 $(LIBDIR)/libcfpu.a
  215. X    -cp fpumath.h $(INCDIR)/fpumath.h
  216. X
  217. clean:
  218. X    rm -f $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ) $(TESTOBJ) core a.out
  219. X
  220. clobber:
  221. X    rm -f $(LIBOBJ) $(LIBFOBJ) $(COMMONOBJ) $(TESTOBJ) libfpu.a \
  222. X    libffpu.a libcfpu.a paranoia fparanoia core a.out mathlib.0[1-6]
  223. X
  224. shar:
  225. X    shar -c -L 55 -a -n mathlib$(VERS) -o mathlib $(ALLSRC) fpumath.h
  226. X
  227. # Dependencies generated using gcc -MM *.c
  228. erf.o : erf.c 
  229. gamma.o : gamma.c 
  230. j0.o : j0.c
  231. j1.o : j1.c
  232. lgamma.o : lgamma.c
  233. X    $(CC) -O -c lgamma.c
  234. ieee_retro.o : ieee_retro.c
  235. nextafter.o : nextafter.c fpumath.h
  236. erff.o : erff.c fpumath.h
  237. gammaf.o : gammaf.c  fpumath.h
  238. j0f.o : j0f.c fpumath.h
  239. j1f.o : j1f.c fpumath.h
  240. nextafterf.o : nextafterf.c fpumath.h
  241. lgammaf.o : lgammaf.c fpumath.h
  242. paranoia.o : paranoia.c 
  243. SHAR_EOF
  244. chmod 0644 Makefile ||
  245. echo 'restore of Makefile failed'
  246. Wc_c="`wc -c < 'Makefile'`"
  247. test 3879 -eq "$Wc_c" ||
  248.     echo 'Makefile: original size 3879, current size' "$Wc_c"
  249. rm -f _shar_wnt_.tmp
  250. fi
  251. # ============= PROBLEMS ==============
  252. if test -f 'PROBLEMS' -a X"$1" != X"-c"; then
  253.     echo 'x - skipping PROBLEMS (File already exists)'
  254.     rm -f _shar_wnt_.tmp
  255. else
  256. > _shar_wnt_.tmp
  257. echo 'x - extracting PROBLEMS (Text)'
  258. sed 's/^X//' << 'SHAR_EOF' > 'PROBLEMS' &&
  259. exp() (and related exp functions) and pow() give rise to LOS exceptions. I
  260. can't think of any other way of calculating these functions. The results I get
  261. agree with the results from a SUN 4 and with Abramowitz and Stegun 
  262. (to 15 places). If anyone knows of another algorithn please let me know.
  263. X
  264. There may still be problems with some functions for pathological arguments.
  265. In some cases I should break the range of the argument up and use
  266. approximations for large and/or small arguments. In fact, the ESIX supplied
  267. exp() is just as bad (good?!?) as mine :-)
  268. SHAR_EOF
  269. chmod 0644 PROBLEMS ||
  270. echo 'restore of PROBLEMS failed'
  271. Wc_c="`wc -c < 'PROBLEMS'`"
  272. test 557 -eq "$Wc_c" ||
  273.     echo 'PROBLEMS: original size 557, current size' "$Wc_c"
  274. rm -f _shar_wnt_.tmp
  275. fi
  276. # ============= README ==============
  277. if test -f 'README' -a X"$1" != X"-c"; then
  278.     echo 'x - skipping README (File already exists)'
  279.     rm -f _shar_wnt_.tmp
  280. else
  281. > _shar_wnt_.tmp
  282. echo 'x - extracting README (Text)'
  283. sed 's/^X//' << 'SHAR_EOF' > 'README' &&
  284. VERSION: This is version 2.0
  285. --------
  286. This is version 2.0 of the alternative 386 math library. It supplants
  287. all previous versions.
  288. X
  289. Version Notes
  290. -------------
  291. The main change between 1.1 and 2.0 is that all functions are now available
  292. in float versions. To use them you need gcc/gas in order to create the
  293. library (libffpu.a) and an ANSI C compiler to compile your source code because
  294. K&R-1 compilers pass all floating point argumnets to functions as doubles and
  295. this is not the desired behaviour.
  296. X
  297. The float function names are formed from the double versions by appending
  298. an 'f', for example sin -> sinf. These functions are fully prototyped in
  299. fpumath.h.
  300. X
  301. Note that jn and yn are not available at this time.
  302. X
  303. Three libraries are now built:
  304. libfpu.a - double version
  305. libffpu.a - float version
  306. libcfpu.a - combined version 
  307. X
  308. The speed difference between floats and doubles is dependent on the
  309. instruction mix - from memory operands or on stack operands. Floats are
  310. faster to load from memory but once on the stack, all operations occur on
  311. the 80 bit temporary real representation.
  312. X
  313. To use the float version effectively you really need an ANSI stdio
  314. implementation.
  315. X
  316. Introduction
  317. ------------
  318. The files in this directory consist of assembler and C source for an 
  319. alternative maths library for Unices (including Xenix) running on an
  320. 80386/80387 combination and using gcc/gas as the compiler system. These
  321. routines are typically from 2 to 10 times faster than those in the 
  322. supplied maths library; assuming that your library, like mine 
  323. (ESIX rev. D and Xenix 2.3.2), does not use the '387 inbuilts to perform 
  324. transcendental calculations.
  325. X
  326. For those of you without a '387 you must use the full emulator and
  327. consequently may not see any speed up. I haven't tried. Under Xenix you
  328. must have a '387---some of the '387 instructions are not emulated. If you
  329. have a '287 your in the same boat; some of the assembler routines won't
  330. work.
  331. X
  332. I have coded the additional IEEE 754 required functions to provide a
  333. conforming double precision implementation. 
  334. X
  335. You require gcc/gas in order to compile the assembler source code.
  336. X
  337. People who are using SUN 386i's (Roadrunner) may also find this useful.
  338. X
  339. I have found that writing good low level numerical code is far harder
  340. than I initially thought it would be. But I've learnt a lot, and I still 
  341. enjoy doing it.
  342. X
  343. File List
  344. ---------
  345. CHANGELOG       atanhf.s        expm1f.s        j0.c            rint.s
  346. COPYING         ceil.s          fabs.s          j0f.c           rintf.s
  347. COPYRIGHT       ceilf.s         fabsf.s         j1.c            scalb.s
  348. Makefile        copysign.s      finite.s        j1f.c           scalbf.s
  349. PROBLEMS        copysignf.s     finitef.s       lgamma.c        setcont.s
  350. README          cos.s           floor.s         lgammaf.c       setinternal.s
  351. TODO            cosf.s          floorf.s        log.s           sin.s
  352. _getsw.s        cosh.s          fmod.s          log10.s         sinf.s
  353. acos.s          coshf.s         fmodf.s         log10f.s        sinh.s
  354. acosf.s         d2dcomb.summ    fpumath.h       log1p.s         sinhf.s
  355. acosh.s         drem.s          gamma.c         log1pf.s        sqrt.s
  356. acoshf.s        dremf.s         gammaf.c        log2.s          sqrtf.s
  357. asin.s          erf.c           hypot.s         log2f.s         sqrtp.s
  358. asinf.s         erff.c          hypotf.s        logb.s          sqrtpf.s
  359. asinh.s         exp.s           ieee_ext.s      logbf.s         tan.s
  360. asinhf.s        exp10.s         ieee_extf.s     logf.s          tanf.s
  361. atan.s          exp10f.s        ieee_retro.c    nextafter.c     tanh.s
  362. atan2.s         exp2.s          ieee_values.s   nextafterf.c    tanhf.s
  363. atan2f.s        exp2f.s         ieee_valuesf.s  paranoia.c
  364. atanf.s         expf.s          infinity.s      pow.s
  365. atanh.s         expm1.s         infinityf.s     powf.s
  366. X
  367. Comments
  368. --------
  369. The additional program paranoia.c attempts to tell how well your floating
  370. point conforms to the IEEE standard. You may like to compare the output
  371. when linked with your existing library and with this one. Paranoia cannot be
  372. compiled using standard 'cc' on Esix 3.2 revision D.
  373. X
  374. The file d2dcomb.summ contains a summary of some test results. See that
  375. file for details.
  376. X
  377. The IEEE specified functions with which you may not be familiar are:
  378. X
  379. double copysign(x,y)
  380. double x,y;
  381. copysign returns x with the sign of y. IEEE denormal will be set if x is a
  382. denormal.
  383. X
  384. double drem(x)
  385. double x;
  386. drem returns the IEEE remainder of x/y - it may be slow.
  387. X
  388. int finite(x)
  389. double x;
  390. finite returns true if -inf < x < inf; false otherwise. Does not raise any
  391. floating point exceptions.
  392. X
  393. double logb(x)
  394. double x;
  395. logb returns the unbiased exponent of its argument.
  396. X
  397. double rint(x)
  398. double x;
  399. rint returns its argument rounded in the prevailing rounding mode.
  400. X
  401. double scalb(x, n)
  402. double x;
  403. int n;
  404. scalb returns x*2^n.
  405. X
  406. Most of the above are just hooks directly into functions provided as single
  407. instructions in the 80387.
  408. X
  409. double infinity()
  410. infinity returns +inf. No exceptions are raised.
  411. X
  412. int isnan(x)
  413. double x;
  414. isnan returns 1 if x is a nan; 0 otherwise. Ieee exceptions are not
  415. affected.
  416. X
  417. int isnormal(x)
  418. double x;
  419. isnormal returns 1 if x is a normalized number; 0 otherwise. Ieee exceptions 
  420. are not affected.
  421. X
  422. int issubnormal(x)
  423. double x;
  424. issubnormal returns 1 if x is not a normalized number; 0 otherwise. Ieee 
  425. exceptions are not affected.
  426. X
  427. int iszero(x)
  428. double x;
  429. iszero returns 1 if x is +/- 0.0; 0 otherwise. Ieee exceptions are not affected.
  430. X
  431. int isinf(x)
  432. double x;
  433. isinf returns 1 if x is +/- inf; 0 otherwise. Ieee exceptions are not affected.
  434. X
  435. int signbit(x)
  436. double x;
  437. signbit return the sign of x. 1 if negative and 0 if positive. Ieee
  438. exceptions are not affected.
  439. X
  440. The is... functions and signbit are in the file ieee_ext.s
  441. X
  442. void ieee_retrospective(f)
  443. FILE *f;
  444. ieee_retrospective prints a list of IEEE exceptions that are currently
  445. active on the 80387 in the file pointed to by f.
  446. X
  447. double
  448. max_normal()
  449. max_normal returns the maximum +ve normalized number. No exceptions are
  450. raised.
  451. X
  452. double
  453. min_normal()
  454. mmin_normal returns the minimum +ve normalized number. No exceptions are
  455. raised.
  456. X
  457. double
  458. max_subnormal()
  459. max_subnormal returns the largest +ve denormalized number. This raises the
  460. denormal exception because of the way floating point numbers are returned.
  461. X
  462. double
  463. min_subnormal()
  464. min_subnormal returns the minimum +ve denormalized number. This raises the
  465. denormal exception because of the way floating point numbers are returned.
  466. X
  467. double
  468. quiet_nan(n)
  469. long n;
  470. quiet_nan returns a quiet nan. The argument is ignored. No exceptions are
  471. raised.
  472. X
  473. double
  474. signaling_nan(n)
  475. long n;
  476. signaling_nan returns a signaling nan. The argument is ignored. This raises
  477. the invalid operation exception because of the way floating point numbers
  478. are returned.
  479. X
  480. double
  481. nextafter(x,y)
  482. double x,y;
  483. nextafter returns the next representable floating point number from x in
  484. the direction of y. Exceptions may be raised depending on the arguments.
  485. X
  486. double erfcx(x)
  487. double x;
  488. erfcx returns x^2 * erfc(x).
  489. X
  490. Additional Functions
  491. --------------------
  492. It is suggested in the book `Programming the 80386' by Crawford and Gelsinger
  493. that all floating point exceptions except `invalid operation' be masked. 
  494. That is the 80387's inbuilt exception handler should be used. If you want this 
  495. behaviour call setinternal() at the start of your code. This function is defined
  496. by:
  497. X
  498. int setinternal()
  499. X
  500. and is declared for your convenience in fpumath.h. The return value is the
  501. current control word. 
  502. X
  503. You can set the control word using setcont(mode). Again, this is defined in
  504. fpumath.h:
  505. X
  506. int setcont(mode)
  507. int mode;
  508. X
  509. This is really only provided to reset the original mode obtained using 
  510. setinternal(). The previous mode is returned.
  511. X
  512. Note that more general forms of these functions are provided in the standard 
  513. library using the fpsetmask and fpgetmask routines.
  514. X
  515. If you use setinternal(), arithmetic operations like 1/0 will return
  516. infinity - inf will be printed if you are printing the output. Note that 0/0 
  517. will still produce a floating point exception; the value of this operation 
  518. is undefined.
  519. X
  520. double sqrtp(x)
  521. double x;
  522. Returns the sqrt of x in 64 bit (double) mode irrespective of the current
  523. precision.
  524. X
  525. Acknowledgements
  526. ----------------
  527. I'd like to thank the developers of the Berkeley Software Distribution who
  528. believe in software freedom and made my job easier by providing C source
  529. for all the functions. I also thank the author of paranoia.c. And also Dan Lau
  530. of Intel. Also thanks go to W.J. Cody of Argonne National Laboratory
  531. (USA).
  532. X
  533. And to everybody who has sent me feedback.
  534. X
  535. X
  536. I have decided to cover the parts I wrote by the GNU General Public Licence 
  537. with my modification (see below), a copy of which is included with this 
  538. distribution.
  539. X
  540. In a nutshell, the Berkeley code is covered by their licence. My code is covered
  541. by GNU's with my modification (see below). This infringes on nobodies rights.
  542. X
  543. Amendment to the GNU General Public Licence (*ONLY* applicable to my code)
  544. -------------------------------------------
  545. I may choose to cover this code by the Free Software Foundation's General
  546. Public Library Licence when it becomes available. At present I make the 
  547. following amendment to the licencing agreement:
  548. X
  549. *******************************************************************************
  550. You may use this library in products which are distributed in binary format
  551. only as long as you provide source code for the library with the distribution
  552. or will supply the source code for the library for a period of 3 years from the
  553. date of providing the binary files.
  554. ********************************************************************************
  555. X
  556. This in no way changes the GNU licence as applicable to other programs.
  557. X
  558. Also I have not modified the GPL it is ditributed in its original form as
  559. it must be.
  560. X
  561. Installation
  562. ------------
  563. Edit the Makefile and change LIBDIR (default /lib) and INCDIR 
  564. (default /usr/include) according to taste.
  565. Type `make' to produce the library (libfpu.a) and paranoia.
  566. Type `make install' to install libfpu.a in $(LIBDIR) and fpumath.h 
  567. in $(INCDIR).
  568. X
  569. X
  570. X
  571. I'm interested in obtaining feedback on the performance of this library and
  572. of being informed of any bugs. I will continue to support this as long as I
  573. have a 32 bit Intel CPU running some form of UNIX and GNU C. My own
  574. development system was ESIX System V.3.2 revision D running on a 20 MHz
  575. 386 (and 387 of course!), gcc 1.37.1 and gas 1.36 with COFF/stab patches.
  576. X
  577. I have also modfied the f2c libF77 to use setcont() and ieee_retrospective()
  578. which makes the whole f2c environment on a 386 look like SUN FORTRAN. If
  579. anyone is interested in these (trivial) changes drop me a line.
  580. X
  581. Please mail comments and bug reports to glenn@qed.physics.su.oz.au.
  582. X
  583. X            Share and Enjoy,
  584. X                Glenn
  585. X
  586. You can also ftp this stuff from suphys.physics.su.oz.au.
  587. X
  588. Glenn Geers                       | "So when it's over, we're back to people.
  589. Department of Theoretical Physics |  Just to prove that human touch can have
  590. The University of Sydney          |  no equal."
  591. Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'
  592. X
  593. Ph: +61 2 692-3241
  594. SHAR_EOF
  595. chmod 0644 README ||
  596. echo 'restore of README failed'
  597. Wc_c="`wc -c < 'README'`"
  598. test 11151 -eq "$Wc_c" ||
  599.     echo 'README: original size 11151, current size' "$Wc_c"
  600. rm -f _shar_wnt_.tmp
  601. fi
  602. # ============= TODO ==============
  603. if test -f 'TODO' -a X"$1" != X"-c"; then
  604.     echo 'x - skipping TODO (File already exists)'
  605.     rm -f _shar_wnt_.tmp
  606. else
  607. > _shar_wnt_.tmp
  608. echo 'x - extracting TODO (Text)'
  609. sed 's/^X//' << 'SHAR_EOF' > 'TODO' &&
  610. 0. Keep modifying the sources to improve performance.
  611. 1. Modify the assembler source so that it works with the standard 'as'.
  612. 2. Add some more esoteric functions (suggestions?).
  613. 3. Provide jn() and yn().
  614. SHAR_EOF
  615. chmod 0644 TODO ||
  616. echo 'restore of TODO failed'
  617. Wc_c="`wc -c < 'TODO'`"
  618. test 204 -eq "$Wc_c" ||
  619.     echo 'TODO: original size 204, current size' "$Wc_c"
  620. rm -f _shar_wnt_.tmp
  621. fi
  622. # ============= d2dcomb.summ ==============
  623. if test -f 'd2dcomb.summ' -a X"$1" != X"-c"; then
  624.     echo 'x - skipping d2dcomb.summ (File already exists)'
  625.     rm -f _shar_wnt_.tmp
  626. else
  627. > _shar_wnt_.tmp
  628. echo 'x - extracting d2dcomb.summ (Text)'
  629. sed 's/^X//' << 'SHAR_EOF' > 'd2dcomb.summ' &&
  630. Comments
  631. --------
  632. These tests were done using routines distributed with the Portable Math 
  633. Library by Fred Fish. I have a sneaking feeling that his exp, atanh and 
  634. atan return bad results (double precision FORTRAN library under TOPS-20) 
  635. since the '387 and SUN 4 results agree. I have partially verified this 
  636. using the extended tables in Abramowitz & Stegun.
  637. X
  638. Results from libfpu.a
  639. dabs:    maximum relative error 0.000000e+00
  640. acos:    maximum relative error 1.395807e-16
  641. acosh:    maximum relative error 1.686042e-16
  642. asin:    maximum relative error 1.378420e-16
  643. asinh:    maximum relative error 1.181272e-13
  644. atan:    maximum relative error 6.406394e-12
  645. atan2:    maximum relative error 2.261728e-15
  646. atanh:    maximum relative error 3.731406e-13
  647. cos:    maximum relative error 0.000000e+00
  648. cosh:    maximum relative error 3.828768e-15
  649. exp:    maximum relative error 4.165239e-13
  650. log:    maximum relative error 0.000000e+00
  651. log10:    maximum relative error 0.000000e+00
  652. sin:    maximum relative error 1.855096e-16
  653. sinh:    maximum relative error 3.828768e-15
  654. sqrt:    maximum relative error 1.252752e-16
  655. tan:    maximum relative error 1.110223e-16
  656. tanh:    maximum relative error 1.201235e-16
  657. X
  658. Results from libfpu.a with extended precision
  659. dabs:    maximum relative error 0.000000e+00
  660. acos:    maximum relative error 1.892807e-16
  661. acosh:    maximum relative error 0.000000e+00
  662. asin:    maximum relative error 1.431811e-16
  663. asinh:    maximum relative error 7.000127e-16
  664. atan:    maximum relative error 6.406394e-12
  665. atan2:    maximum relative error 2.261728e-15
  666. atanh:    maximum relative error 3.731406e-13
  667. cos:    maximum relative error 0.000000e+00
  668. cosh:    maximum relative error 0.000000e+00
  669. exp:    maximum relative error 4.176205e-13
  670. log:    maximum relative error 0.000000e+00
  671. log10:    maximum relative error 0.000000e+00
  672. sin:    maximum relative error 1.855096e-16
  673. sinh:    maximum relative error 0.000000e+00
  674. sqrt:    maximum relative error 1.252752e-16
  675. tan:    maximum relative error 1.110223e-16
  676. tanh:    maximum relative error 1.226565e-16
  677. X
  678. Results form libm.a (Esix standard) - deleted functions do not exist
  679. dabs:    maximum relative error 0.000000e+00
  680. acos:    maximum relative error 1.892807e-16
  681. asin:    maximum relative error 1.431811e-16
  682. atan:    maximum relative error 6.406394e-12
  683. atan2:    maximum relative error 2.261728e-15
  684. cos:    maximum relative error 0.000000e+00
  685. cosh:    maximum relative error 3.828768e-15
  686. exp:    maximum relative error 4.165239e-13
  687. log:    maximum relative error 0.000000e+00
  688. log10:    maximum relative error 0.000000e+00
  689. sin:    maximum relative error 1.855096e-16
  690. sinh:    maximum relative error 3.828768e-15
  691. sqrt:    maximum relative error 1.252752e-16
  692. tan:    maximum relative error 1.110223e-16
  693. tanh:    maximum relative error 1.201235e-16
  694. X
  695. Results form libm.a (SUN 4)
  696. dabs:    maximum relative error 0.000000e+00
  697. acos:    maximum relative error 0.000000e+00
  698. acosh:    maximum relative error 1.686042e-16
  699. asin:    maximum relative error 1.378420e-16
  700. asinh:    maximum relative error 7.000127e-16
  701. atan:    maximum relative error 6.406394e-12
  702. atan2:    maximum relative error 2.261728e-15
  703. atanh:    maximum relative error 3.731406e-13
  704. cos:    maximum relative error 1.156277e-16
  705. cosh:    maximum relative error 1.651640e-16
  706. exp:    maximum relative error 4.176205e-13
  707. log:    maximum relative error 0.000000e+00
  708. log10:    maximum relative error 0.000000e+00
  709. sin:    maximum relative error 1.855096e-16
  710. sinh:    maximum relative error 0.000000e+00
  711. sqrt:    maximum relative error 1.252753e-16
  712. tan:    maximum relative error 1.425732e-16
  713. tanh:    maximum relative error 0.000000e+00
  714. SHAR_EOF
  715. chmod 0644 d2dcomb.summ ||
  716. echo 'restore of d2dcomb.summ failed'
  717. Wc_c="`wc -c < 'd2dcomb.summ'`"
  718. test 3424 -eq "$Wc_c" ||
  719.     echo 'd2dcomb.summ: original size 3424, current size' "$Wc_c"
  720. rm -f _shar_wnt_.tmp
  721. fi
  722. # ============= fpumath.h ==============
  723. if test -f 'fpumath.h' -a X"$1" != X"-c"; then
  724.     echo 'x - skipping fpumath.h (File already exists)'
  725.     rm -f _shar_wnt_.tmp
  726. else
  727. > _shar_wnt_.tmp
  728. echo 'x - extracting fpumath.h (Text)'
  729. sed 's/^X//' << 'SHAR_EOF' > 'fpumath.h' &&
  730. /*
  731. ** This file is covered by the GNU General Public Licence
  732. ** 
  733. ** This file should be installed somewhere in your standard include
  734. ** search path and given the name fpumath.h.
  735. */
  736. X
  737. /*
  738. ** This is a good way of flagging whether 
  739. ** you are using the alternate stuff
  740. */
  741. #ifndef    __FPUMATH_H
  742. #define    __FPUMATH_H
  743. X
  744. #define __FPU__
  745. X
  746. X
  747. #if defined(__GNUC__) || defined(__STDC__)
  748. X
  749. extern double j0(double), j1(double), jn(int,double), y0(double), y1(double);
  750. extern double yn(int,double);
  751. extern double erf(double), erfc(double), erfcx(double);
  752. extern double exp(double), log(double), log10(double), pow(double,double);
  753. extern double sqrt(double);
  754. extern double expm1(double), log1p(double), exp10(double);
  755. extern double exp2(double), log2(double);
  756. extern double floor(double), ceil(double), fabs(double);
  757. extern double lgamma(double), gamma(double);
  758. extern double sinh(double), cosh(double), tanh(double);
  759. extern double asinh(double), acosh(double), atanh(double);
  760. extern double sin(double), cos(double), tan(double), asin(double);
  761. extern double acos(double), atan(double), atan2(double,double);
  762. extern double hypot(double,double);
  763. extern double sqrtp(double);
  764. X
  765. #else
  766. X
  767. extern double j0(), j1(), jn(), y0(), y1(), yn();
  768. extern double erf(), erfc(), erfcx();
  769. extern double exp(), log(), log10(), pow(), sqrt();
  770. extern double expm1(), log1p(), exp10();
  771. extern double exp2(), log2();
  772. extern double floor(), ceil(), fabs();
  773. extern double lgamma(), gamma();
  774. extern double sinh(), cosh(), tanh();
  775. extern double asinh(), acosh(), atanh();
  776. extern double sin(), cos(), tan(), asin();
  777. extern double acos(), atan(), atan2(), hypot();
  778. extern double sqrtp();
  779. X
  780. #endif
  781. X
  782. /*
  783. ** IEEE functions
  784. */
  785. X
  786. #if defined(__GNUC__) || defined(__STDC__)
  787. X
  788. extern double rint(double), drem(double,double), scalb(double,int);
  789. extern double logb(double), copysign(double,double);
  790. extern double infinity(void), nextafter(double,double);
  791. extern int isnan(double), iszero(double), isinf(double), isnormal(double);
  792. extern int issubnormal(double);
  793. extern int signbit(double), finite(double);
  794. extern double max_normal(void), min_normal(void), min_subnormal(void);
  795. extern double max_subnormal(void);
  796. extern double quiet_nan(long int),signaling_nan(long int);
  797. X
  798. #else
  799. X
  800. extern double rint(), drem(), scalb(), logb(), copysign();
  801. extern double infinity(), nextafter();
  802. extern int isnan(), iszero(), isinf(), isnormal(), issubnormal();
  803. extern int signbit(), finite();
  804. extern double max_normal(), min_normal(), min_subnormal(), max_subnormal();
  805. extern double quiet_nan(),signaling_nan();
  806. X
  807. #endif
  808. X
  809. /*
  810. ** Float versions
  811. */
  812. #if defined(__GNUC__) || defined(__STDC__)
  813. X
  814. extern float j0f(float), j1f(float), jnf(int,float), y0f(float);
  815. extern float y1f(float), ynf(int,float);
  816. extern float erff(float), erfcf(float), erfcxf(float), sqrtf(float);
  817. extern float expf(float), logf(float), log10f(float), powf(float, float);
  818. extern float expm1f(float), log1pf(float), exp10f(float);
  819. extern float exp2f(float), log2f(float);
  820. extern float floorf(float), ceilf(float), fabsf(float);
  821. extern float lgammaf(float), gammaf(float);
  822. extern float sinhf(float), coshf(float), tanhf(float);
  823. extern float asinhf(float), acoshf(float), atanhf(float);
  824. extern float sinf(float), cosf(float), tanf(float), asinf(float);
  825. extern float acosf(float), atanf(float), atan2f(float, float);
  826. extern float hypotf(float, float);
  827. extern float sqrtpf(float);
  828. X
  829. extern float rintf(float), dremf(float), scalbf(float, int);
  830. extern float logbf(float), copysignf(float, float);
  831. extern float infinityf(void), nextafterf(float, float);
  832. extern int isnanf(float), iszerof(float), isinff(float);
  833. extern int isnormalf(float), issubnormalf(float);
  834. extern int signbitf(float), finitef(float);
  835. extern float max_normalf(void), min_normalf(void);
  836. extern float min_subnormalf(void), max_subnormalf(void);
  837. extern float quiet_nanf(long int), signaling_nanf(long int);
  838. X
  839. #endif
  840. X
  841. /*
  842. ** Extensions
  843. */
  844. X
  845. #if defined(__GNUC__) || defined(__STDC__)
  846. X
  847. extern int setinternal(void), setcont(int);
  848. X
  849. #else
  850. X
  851. extern int setinternal(), setcont();
  852. X
  853. #endif
  854. X
  855. /*
  856. ** The argument of ieee_retrospective is a FILE * 
  857. ** but we don't want to include stdio.h.
  858. */
  859. X
  860. extern void ieee_retrospective();
  861. X
  862. /*
  863. ** Useful defines for setcont
  864. */
  865. X
  866. #define MASK_ALL0 0x107f    /* single precision */
  867. #define MASK_ALL 0x127f        /* double precision */
  868. #define MASK_ALL1 0x137f     /* extra precision */
  869. X
  870. #define __infinity infinity
  871. X
  872. #define M_E    2.7182818284590452354
  873. #define M_LOG2E    1.4426950408889634074
  874. #define M_LOG10E    0.43429448190325182765
  875. #define M_LN2    0.69314718055994530942
  876. #define M_LN10    2.30258509299404568402
  877. #define M_PI    3.14159265358979323846
  878. #define M_PI_2    1.57079632679489661923
  879. #define M_PI_4    0.78539816339744830962
  880. #define M_1_PI    0.31830988618379067154
  881. #define M_2_PI    0.63661977236758134308
  882. #define M_2_SQRTPI    1.12837916709551257390
  883. #define M_SQRT2    1.41421356237309504880
  884. #define M_SQRT1_2    0.70710678118654752440
  885. X
  886. #ifndef IEEE
  887. #define MAXFLOAT    ((float)3.40282346638528860e+38)
  888. #define    MINFLOAT    ((float)1.40129846432481707e-45)
  889. #else
  890. #define    MAXFLOAT    (max_normalf())
  891. #define    MINFLOAT    (min_subnormalf())
  892. #endif
  893. X
  894. #ifndef IEEE
  895. #define    MAXDOUBLE    1.79769313486231570e+308  /* wrong in math.h */
  896. #define    MINDOUBLE    4.94065645841246544e-324
  897. #else
  898. #define    MAXDOUBLE    (max_normal())
  899. #define    MINDOUBLE    (min_subnormal())
  900. #endif
  901. X
  902. #define ULP        1.1102230246251565e-16
  903. #define ULPF        5.9604644775390625e-08
  904. X
  905. #define    HUGE_VAL    MAXFLOAT
  906. X
  907. #ifdef IEEE
  908. #define HUGE    (__infinity())
  909. #else
  910. #define HUGE     MAXDOUBLE
  911. #endif
  912. X
  913. #endif
  914. SHAR_EOF
  915. chmod 0644 fpumath.h ||
  916. echo 'restore of fpumath.h failed'
  917. Wc_c="`wc -c < 'fpumath.h'`"
  918. test 5452 -eq "$Wc_c" ||
  919.     echo 'fpumath.h: original size 5452, current size' "$Wc_c"
  920. rm -f _shar_wnt_.tmp
  921. fi
  922. rm -f _shar_seq_.tmp
  923. echo You have unpacked the last part
  924. exit 0
  925. --
  926. Glenn Geers                       | "So when it's over, we're back to people.
  927. Department of Theoretical Physics |  Just to prove that human touch can have
  928. The University of Sydney          |  no equal."
  929. Sydney NSW 2006 Australia         |  - Basia Trzetrzelewska, 'Prime Time TV'
  930.