add JavaScript Interpreter

This commit is contained in:
geniusgogo
2014-02-25 01:47:49 +08:00
parent 7255137b0a
commit 121bb5fcdf
115 changed files with 29947 additions and 13 deletions

View File

@@ -0,0 +1,46 @@
Contributing
===========
Thanks for thinking about contributing to Espruino! Anything you can add is hugely appreciated, but please can you follow a few simple rules:
* Keep the same coding style (See **Coding Style** below)
* Ensure that you are not contributing someone else's code, and that you are willing to add your code under Espruino's MPL Licence
* Make sure that what you do doesn't break the Espruino board or the other boards we build for. We can't check all the boards for every commit, so if you break something you'll annoy a whole bunch of people.
* Be aware that Espruino is designed for Microcontrollers - with very low amounts of flash and memory. Both are at a premium so don't statically allocate variables or do other stuff that will use up memory.
* Avoid randomly adding newlines, spaces, refactoring everything or renaming things to your own personal style (some things really could do with renaming, but please check first or we may reject your pull request)
* Don't add a whole bunch of indirection/abstraction for the sake of it - it'll probably just use of more of our precious memory.
* If you add a new API, try and make it familiar to Arduino/JavaScript users.
Target Areas
-----------
We'll keep the outstanding issues in GitHub's issue list, but general stuff that would really help us is:
* Tests. If something doesn't work, please make a test for it. Even if you don't fix it it'll help others greatly. Bonus points if it's in a pull request :)
* Documentation. Improving the documentation (either the EspruinoDocs project, or the auto-generated reference) would be fantastic.
* Duplication. If the same code is used for multiple platforms, try and make sure it's shared, not duplicated.
* Remove hard-coded stuff. Various things like the SPI filesystem are still hard-coded with ifdefs for each board - we want all that stuff to be generated from build_platform_info.py
* Speed. There are a few areas this could be improved - but please benchmark what you're doing both before and afterwards on the Espruino board to check that what you've done helps
* Memory Usage. Both RAM and Flash are at a premium. Ways of reducing this (including stack usage) and making usage more efficient are really appreciated.
* JavaScript compliance (without affecting speed or memory usage too much).
Contributing
-----------
* Please RUN THE TESTS to check that there are no regressions
* Issue us a pull request to [www.github.com/espruino] via GitHub
* Please keep each request small (just include one fix per request)
Coding Style
-----------
The rough coding style is as follows, but you should get a good idea from the code. If we've missed anything obvious please let us know!
* 2 Spaces for indents
* Open curly braces on the same line
* No Tabs used
* Use bool for booleans - not int
* ```//``` comments for single lines, ```/* ... */``` for multiple lines
* Half-hearted Doxygen compatibility: use ```///<``` for function declaration documentation (if on same line), and ```/** ... */``` if doing it right before a function
* Use new lines sparingly (only where it really makes sense)

515
components/external/espruino/ChangeLog vendored Normal file

File diff suppressed because it is too large Load Diff

378
components/external/espruino/LICENSE vendored Normal file
View File

@@ -0,0 +1,378 @@
All files in this package are Copyright 2013 Gordon Williams, Pur3 Ltd unless
otherwise noted.
-------------------------------------------------------------------------------
Mozilla Public License Version 2.0
==================================
1. Definitions
--------------
1.1. "Contributor"
means each individual or legal entity that creates, contributes to
the creation of, or owns Covered Software.
1.2. "Contributor Version"
means the combination of the Contributions of others (if any) used
by a Contributor and that particular Contributor's Contribution.
1.3. "Contribution"
means Covered Software of a particular Contributor.
1.4. "Covered Software"
means Source Code Form to which the initial Contributor has attached
the notice in Exhibit A, the Executable Form of such Source Code
Form, and Modifications of such Source Code Form, in each case
including portions thereof.
1.5. "Incompatible With Secondary Licenses"
means
(a) that the initial Contributor has attached the notice described
in Exhibit B to the Covered Software; or
(b) that the Covered Software was made available under the terms of
version 1.1 or earlier of the License, but not also under the
terms of a Secondary License.
1.6. "Executable Form"
means any form of the work other than Source Code Form.
1.7. "Larger Work"
means a work that combines Covered Software with other material, in
a separate file or files, that is not Covered Software.
1.8. "License"
means this document.
1.9. "Licensable"
means having the right to grant, to the maximum extent possible,
whether at the time of the initial grant or subsequently, any and
all of the rights conveyed by this License.
1.10. "Modifications"
means any of the following:
(a) any file in Source Code Form that results from an addition to,
deletion from, or modification of the contents of Covered
Software; or
(b) any new file in Source Code Form that contains any Covered
Software.
1.11. "Patent Claims" of a Contributor
means any patent claim(s), including without limitation, method,
process, and apparatus claims, in any patent Licensable by such
Contributor that would be infringed, but for the grant of the
License, by the making, using, selling, offering for sale, having
made, import, or transfer of either its Contributions or its
Contributor Version.
1.12. "Secondary License"
means either the GNU General Public License, Version 2.0, the GNU
Lesser General Public License, Version 2.1, the GNU Affero General
Public License, Version 3.0, or any later versions of those
licenses.
1.13. "Source Code Form"
means the form of the work preferred for making modifications.
1.14. "You" (or "Your")
means an individual or a legal entity exercising rights under this
License. For legal entities, "You" includes any entity that
controls, is controlled by, or is under common control with You. For
purposes of this definition, "control" means (a) the power, direct
or indirect, to cause the direction or management of such entity,
whether by contract or otherwise, or (b) ownership of more than
fifty percent (50%) of the outstanding shares or beneficial
ownership of such entity.
2. License Grants and Conditions
--------------------------------
2.1. Grants
Each Contributor hereby grants You a world-wide, royalty-free,
non-exclusive license:
(a) under intellectual property rights (other than patent or trademark)
Licensable by such Contributor to use, reproduce, make available,
modify, display, perform, distribute, and otherwise exploit its
Contributions, either on an unmodified basis, with Modifications, or
as part of a Larger Work; and
(b) under Patent Claims of such Contributor to make, use, sell, offer
for sale, have made, import, and otherwise transfer either its
Contributions or its Contributor Version.
2.2. Effective Date
The licenses granted in Section 2.1 with respect to any Contribution
become effective for each Contribution on the date the Contributor first
distributes such Contribution.
2.3. Limitations on Grant Scope
The licenses granted in this Section 2 are the only rights granted under
this License. No additional rights or licenses will be implied from the
distribution or licensing of Covered Software under this License.
Notwithstanding Section 2.1(b) above, no patent license is granted by a
Contributor:
(a) for any code that a Contributor has removed from Covered Software;
or
(b) for infringements caused by: (i) Your and any other third party's
modifications of Covered Software, or (ii) the combination of its
Contributions with other software (except as part of its Contributor
Version); or
(c) under Patent Claims infringed by Covered Software in the absence of
its Contributions.
This License does not grant any rights in the trademarks, service marks,
or logos of any Contributor (except as may be necessary to comply with
the notice requirements in Section 3.4).
2.4. Subsequent Licenses
No Contributor makes additional grants as a result of Your choice to
distribute the Covered Software under a subsequent version of this
License (see Section 10.2) or under the terms of a Secondary License (if
permitted under the terms of Section 3.3).
2.5. Representation
Each Contributor represents that the Contributor believes its
Contributions are its original creation(s) or it has sufficient rights
to grant the rights to its Contributions conveyed by this License.
2.6. Fair Use
This License is not intended to limit any rights You have under
applicable copyright doctrines of fair use, fair dealing, or other
equivalents.
2.7. Conditions
Sections 3.1, 3.2, 3.3, and 3.4 are conditions of the licenses granted
in Section 2.1.
3. Responsibilities
-------------------
3.1. Distribution of Source Form
All distribution of Covered Software in Source Code Form, including any
Modifications that You create or to which You contribute, must be under
the terms of this License. You must inform recipients that the Source
Code Form of the Covered Software is governed by the terms of this
License, and how they can obtain a copy of this License. You may not
attempt to alter or restrict the recipients' rights in the Source Code
Form.
3.2. Distribution of Executable Form
If You distribute Covered Software in Executable Form then:
(a) such Covered Software must also be made available in Source Code
Form, as described in Section 3.1, and You must inform recipients of
the Executable Form how they can obtain a copy of such Source Code
Form by reasonable means in a timely manner, at a charge no more
than the cost of distribution to the recipient; and
(b) You may distribute such Executable Form under the terms of this
License, or sublicense it under different terms, provided that the
license for the Executable Form does not attempt to limit or alter
the recipients' rights in the Source Code Form under this License.
3.3. Distribution of a Larger Work
You may create and distribute a Larger Work under terms of Your choice,
provided that You also comply with the requirements of this License for
the Covered Software. If the Larger Work is a combination of Covered
Software with a work governed by one or more Secondary Licenses, and the
Covered Software is not Incompatible With Secondary Licenses, this
License permits You to additionally distribute such Covered Software
under the terms of such Secondary License(s), so that the recipient of
the Larger Work may, at their option, further distribute the Covered
Software under the terms of either this License or such Secondary
License(s).
3.4. Notices
You may not remove or alter the substance of any license notices
(including copyright notices, patent notices, disclaimers of warranty,
or limitations of liability) contained within the Source Code Form of
the Covered Software, except that You may alter any license notices to
the extent required to remedy known factual inaccuracies.
3.5. Application of Additional Terms
You may choose to offer, and to charge a fee for, warranty, support,
indemnity or liability obligations to one or more recipients of Covered
Software. However, You may do so only on Your own behalf, and not on
behalf of any Contributor. You must make it absolutely clear that any
such warranty, support, indemnity, or liability obligation is offered by
You alone, and You hereby agree to indemnify every Contributor for any
liability incurred by such Contributor as a result of warranty, support,
indemnity or liability terms You offer. You may include additional
disclaimers of warranty and limitations of liability specific to any
jurisdiction.
4. Inability to Comply Due to Statute or Regulation
---------------------------------------------------
If it is impossible for You to comply with any of the terms of this
License with respect to some or all of the Covered Software due to
statute, judicial order, or regulation then You must: (a) comply with
the terms of this License to the maximum extent possible; and (b)
describe the limitations and the code they affect. Such description must
be placed in a text file included with all distributions of the Covered
Software under this License. Except to the extent prohibited by statute
or regulation, such description must be sufficiently detailed for a
recipient of ordinary skill to be able to understand it.
5. Termination
--------------
5.1. The rights granted under this License will terminate automatically
if You fail to comply with any of its terms. However, if You become
compliant, then the rights granted under this License from a particular
Contributor are reinstated (a) provisionally, unless and until such
Contributor explicitly and finally terminates Your grants, and (b) on an
ongoing basis, if such Contributor fails to notify You of the
non-compliance by some reasonable means prior to 60 days after You have
come back into compliance. Moreover, Your grants from a particular
Contributor are reinstated on an ongoing basis if such Contributor
notifies You of the non-compliance by some reasonable means, this is the
first time You have received notice of non-compliance with this License
from such Contributor, and You become compliant prior to 30 days after
Your receipt of the notice.
5.2. If You initiate litigation against any entity by asserting a patent
infringement claim (excluding declaratory judgment actions,
counter-claims, and cross-claims) alleging that a Contributor Version
directly or indirectly infringes any patent, then the rights granted to
You by any and all Contributors for the Covered Software under Section
2.1 of this License shall terminate.
5.3. In the event of termination under Sections 5.1 or 5.2 above, all
end user license agreements (excluding distributors and resellers) which
have been validly granted by You or Your distributors under this License
prior to termination shall survive termination.
************************************************************************
* *
* 6. Disclaimer of Warranty *
* ------------------------- *
* *
* Covered Software is provided under this License on an "as is" *
* basis, without warranty of any kind, either expressed, implied, or *
* statutory, including, without limitation, warranties that the *
* Covered Software is free of defects, merchantable, fit for a *
* particular purpose or non-infringing. The entire risk as to the *
* quality and performance of the Covered Software is with You. *
* Should any Covered Software prove defective in any respect, You *
* (not any Contributor) assume the cost of any necessary servicing, *
* repair, or correction. This disclaimer of warranty constitutes an *
* essential part of this License. No use of any Covered Software is *
* authorized under this License except under this disclaimer. *
* *
************************************************************************
************************************************************************
* *
* 7. Limitation of Liability *
* -------------------------- *
* *
* Under no circumstances and under no legal theory, whether tort *
* (including negligence), contract, or otherwise, shall any *
* Contributor, or anyone who distributes Covered Software as *
* permitted above, be liable to You for any direct, indirect, *
* special, incidental, or consequential damages of any character *
* including, without limitation, damages for lost profits, loss of *
* goodwill, work stoppage, computer failure or malfunction, or any *
* and all other commercial damages or losses, even if such party *
* shall have been informed of the possibility of such damages. This *
* limitation of liability shall not apply to liability for death or *
* personal injury resulting from such party's negligence to the *
* extent applicable law prohibits such limitation. Some *
* jurisdictions do not allow the exclusion or limitation of *
* incidental or consequential damages, so this exclusion and *
* limitation may not apply to You. *
* *
************************************************************************
8. Litigation
-------------
Any litigation relating to this License may be brought only in the
courts of a jurisdiction where the defendant maintains its principal
place of business and such litigation shall be governed by laws of that
jurisdiction, without reference to its conflict-of-law provisions.
Nothing in this Section shall prevent a party's ability to bring
cross-claims or counter-claims.
9. Miscellaneous
----------------
This License represents the complete agreement concerning the subject
matter hereof. If any provision of this License is held to be
unenforceable, such provision shall be reformed only to the extent
necessary to make it enforceable. Any law or regulation which provides
that the language of a contract shall be construed against the drafter
shall not be used to construe this License against a Contributor.
10. Versions of the License
---------------------------
10.1. New Versions
Mozilla Foundation is the license steward. Except as provided in Section
10.3, no one other than the license steward has the right to modify or
publish new versions of this License. Each version will be given a
distinguishing version number.
10.2. Effect of New Versions
You may distribute the Covered Software under the terms of the version
of the License under which You originally received the Covered Software,
or under the terms of any subsequent version published by the license
steward.
10.3. Modified Versions
If you create software not governed by this License, and you want to
create a new license for such software, you may create and use a
modified version of this License if you rename the license and remove
any references to the name of the license steward (except to note that
such modified license differs from this License).
10.4. Distributing Source Code Form that is Incompatible With Secondary
Licenses
If You choose to distribute Source Code Form that is Incompatible With
Secondary Licenses under the terms of this version of the License, the
notice described in Exhibit B of this License must be attached.
Exhibit A - Source Code Form License Notice
-------------------------------------------
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
If it is not possible or desirable to put the notice in a particular
file, then You may include the notice in a location (such as a LICENSE
file in a relevant directory) where a recipient would be likely to look
for such a notice.
You may add additional accurate notices of copyright ownership.
Exhibit B - "Incompatible With Secondary Licenses" Notice
---------------------------------------------------------
This Source Code Form is "Incompatible With Secondary Licenses", as
defined by the Mozilla Public License, v. 2.0.

969
components/external/espruino/Makefile vendored Normal file

File diff suppressed because it is too large Load Diff

161
components/external/espruino/README.md vendored Normal file
View File

@@ -0,0 +1,161 @@
Espruino JavaScript for Microcontrollers
========================================
<pre>
_____ _
| __|___ ___ ___ _ _|_|___ ___
| __|_ -| . | _| | | | | . |
|_____|___| _|_| |___|_|_|_|___|
|_|
</pre>
http://www.espruino.com
**NOTE:** This software is beta and is provided as-is, and won't be considered even remotely final until we've released the Espruino Board. As such, don't expect support, and do expect it to change rapidly and without warning. Build your own documentation (see **Building**), as the API may be different from the one described on the Espruino website.
The KickStarter campaign said the Espruino Board will have some things which this repository does not yet have (like working CC3000 support). These are works in progress and should be done by the time you get your board (or will be available as a software update).
About
-----
It'd probably help to read the [FAQ](http://www.espruino.com/FAQ), and specifically the page about [Performance](http://www.espruino.com/Performance) as it contains information about how Espruino itself works.
There's also the auto-generated [Reference](http://www.espruino.com/Reference) for JavaScript commands as well as the [Tutorials](http://www.espruino.com/Tutorials) on the website. However please note that this repository is under heavy development, and the documentation on the Espruino website will match the version [available for download](http://www.espruino.com/Download) but **not** the latest version from Git.
License
-------
Please see the [LICENSE](LICENSE) file
Found a Bug?
------------
Please check that:
* It hasn't [already been found](https://github.com/espruino/Espruino/issues) or [been covered on our forum](www.espruino.com/Forum)
* You're not just looking at outdated documentation (See the [Building](#Building) section to see how to build documentation)
Please [submit bugs](https://github.com/espruino/Espruino/issues) with clear steps to reproduce them (and ideally a test case for the ```tests``` directory), and if at all possible try and include a patch to fix them. Please be aware that we have a whole bunch of outstanding issues (some quite large), so if you report something (especially if it doesn't contain a test or a pull request) it may not be fixed for quite some time.
Contributing
------------
Please see [CONTRIBUTING.md](CONTRIBUTING.md)
Current State
-------------
You can download binaries from http://www.espruino.com/Download (these aren't the latest, but are more likely to work with your board)
Please note that this is BETA. We've been working hard on the Espruino Board support but we haven't had time to check the other boards properly.
* Espruino Board - working
* Linux - working
* STM32VLDISCOVERY - WORKING
* STM32F3DISCOVERY - WORKING
* STM32F4DISCOVERY - WORKING
* STM32F429IDISCOVERY - not working, currently no LCD support
* HY STM32 2.4" - NOT WORKING - appears to crash after startup
* HY STM32 2.8" - WORKING, but screen is not black at startup
* HY STM32 3.2" - WORKING
* Olimexino - WORKING
* Carambola - ?
* Raspberry Pi - WORKING
* Sony SmartWatch - USB VCP support still needed
* MBed platforms - have not worked for a while - hardware wrapper still needed
* Arduino - has never worked. Compiles but doesn't even get past init
* LC-TECH STM32F103RBT6 - WORKING, but with some issues (LED inverted logic, BTN needs pullup to work)
Using
-----
If you're using Espruino for your own personal projects - go ahead, we hope you have fun - and please let us know what you do with it on http://www.espruino.com/Forum!
However if you're planning on selling the Espruino software on your own board, please talk to us:
* Read the terms of the MPLv2 Licence that Espruino is distributed under, and make sure you comply with it
* You won't be able to call your board 'Espruino' but you must explain clearly that it uses 'Espruino' internally (we own the trademark)
* If you're profiting from Espruino without contributing anything back, we won't support you (or your users)
Building
--------
Espruino is easy to build under Linux, and it is possible to build under MacOS. We'd strongly suggest that you DO NOT TRY AND BUILD UNDER WINDOWS, and instead use a Virtual Machine. There's a good post on this here: http://forum.espruino.com/conversations/151
We suggest that you use the CodeSourcery GCC compiler, but paths in Makefile may need changing...
``` BOARDNAME=1 RELEASE=1 make```
* See the top of Makefile for board names
* Without `RELEASE=1`, assertions are kept in the code (which is good for debugging, bad for performance + code size)
* `BOARDNAME=1 RELEASE=1 make serialflash` will flash to /dev/ttyUSB0 using the STM32 serial bootloader (what's needed for Espruino + HY boards)
* `BOARDNAME=1 RELEASE=1 make flash` will flash using st-flash if discovery, or maple bootloader if using that board
You can build documentation by running:
``` python scripts/build_docs.py ```
This will create a file called ```functions.html```
Directories and Files
---------------------
* `ChangeLog`: What's new
* `TODO`: List of things to do
* `boards/`: Information on boards, used to auto-generate a lot of the code
* `code/`: Example JavaScript code
* `gen/`: Auto-Generated Source Files
* `libs/`: Optional libraries to include in Espruino (Math, Filesystem, Graphics, etc)
* `linker/`: Linker files for various processors
* `misc/`: random other stuff
* `scripts/`: Scripts for generating files in gen, and for analysing code/compilation/etc
* `src/`: Main source code
* `targetlibs/`: Libraries for targeted architectures
* `targets/`: Specific code for targeted architectures
* `tests/`: JavaScript Testcases
* `benchmark/`: JavaScript Benchmarks
* `dist_*`: files to be copied into distribution zip file
Adding more devices
-------------------
Currently there are a bunch of different files to modify. Eventually the plan is to fit everything into boards/BOARDNAME.py and to auto-generate the rest of the config files.
* Most build options handled in `Makefile`
* Extra libraries like USB/LCD/filesystem in `Makefile`
* Linker Scripts are in `linker/`
* `boards/*.py` files handle loading the list of available pins so the relevant headers + docs can be created
* Processor-specific code in `targets/stm32`, `targets/linux`, etc.
* Processor-specific libs in `targetlibs/foo`
* `src/jshardware.h` is effectively a simple abstraction layer for SPI/I2C/etc
* `targets/stm32/jshardware.c` also has flash-size-specific defines
* `libs/fat_sd` and `libs/lcd` still have some device-specific defines in too
Adding libraries
-------------------
* Create `jswrap_mylib.c/h` in `libs/`
* Create library functions (see examples in other jswrap files, also the comments in `scripts/common.py`)
Arduino Compile (beta)
----------------------
* Ensure that `targets/arduino/utility` is symlinked to `src`
* Symlink `...arduino_workspace/libraries/Espruino` to `targets/arduino`
Cross Compile for Raspberry Pi
------------------------------
```
cd targetlibs
mkdir raspberrypi
cd raspberrypi
git clone git://github.com/raspberrypi/tools.git
sudo apt-get install ia32-libs
```
Cross Compile for Carambola (OpenWRT)
-------------------------------------
* Follow instructions at <https://github.com/8devices/carambola> to set toolchain up in ```~/workspace/carambola```
* Run ```CARAMBOLA=1 make```

View File

@@ -0,0 +1,2 @@
这是从Espruino移植过来的javascript解释器。
这个组件依赖newlib。需要在rtconfig.h中开启RT_USING_JS宏来开启此组件。

View File

@@ -1,19 +1,58 @@
import rtconfig
Import('RTT_ROOT')
from building import *
import os
cwd = GetCurrentDir()
# source files
src = Split('''
''')
path = [RTT_ROOT + r'/components/external/Espruino']
path = path + [RTT_ROOT + r'/components/external/Espruino/src/']
path = path + [RTT_ROOT + r'/components/external/Espruino/gen/']
path = path + [RTT_ROOT + r'/components/external/Espruino/libs/']
path = path + [RTT_ROOT + r'/components/external/Espruino/targets/rtthread/']
ESPRUINO_SRC_PATH = cwd + '/src'
if GetDepend('RT_USING_JS') and not os.path.exists(ESPRUINO_SRC_PATH):
print '================ERROR============================'
print 'Please get espruino source files and put them under espruino folder'
print '================================================='
exit(0)
src = Split("""
src/jsdevices.c
src/jslex.c
src/jswrap_array.c
src/jswrap_json.c
src/jswrap_process.c
src/jsutils.c
src/jswrap_interactive.c
src/jswrap_onewire.c
src/jswrap_string.c
src/jsparse.c
src/jswrap_arraybuffer.c
src/jswrap_modules.c
src/jswrap_serial.c
src/jsinteractive.c
src/jsvar.c
src/jswrap_io.c
src/jswrap_pin.c
src/jspin.c
src/jswrap_functions.c
src/jswrap_object.c
src/jswrap_spi_i2c.c
"""
)
CPPPATH = [cwd + '/src', cwd + '/gen', cwd + '/libs', cwd + '/targets/rtthread']
group = DefineGroup('Espruino', src, depend = ['RT_USING_JS'], CPPPATH = CPPPATH)
gen = Split("""
gen/jspininfo.c
gen/jswrapper.c
"""
)
target = Split("""
targets/rtthread/espruino.c
targets/rtthread/jshardware.c
"""
)
libs = Split(
"""
libs/jswrap_math.c
"""
)
src += gen + target + libs
group = DefineGroup('Espruino', src, depend = ['RT_USING_JS', 'RT_USING_NEWLIB'], CPPPATH = path)
Return('group')

View File

@@ -0,0 +1,109 @@
Espruino uses a few libraries of code that have been generously given away for
free. Their licences are included below.
-------------------------------------------------------------------------
FATFS
-------------------------------------------------------------------------
/*---------------------------------------------------------------------------/
/ FatFs - FAT file system module include file R0.07c (C)ChaN, 2009
/----------------------------------------------------------------------------/
/ FatFs module is an open source software to implement FAT file system to
/ small embedded systems. This is a free software and is opened for education,
/ research and commercial developments under license policy of following terms.
/
/ Copyright (C) 2009, ChaN, all right reserved.
/
/ * The FatFs module is a free software and there is NO WARRANTY.
/ * No restriction on use. You can use, modify and redistribute it for
/ personal, non-profit or commercial product UNDER YOUR RESPONSIBILITY.
/ * Redistributions of source code must retain the above copyright notice.*/
-------------------------------------------------------------------------
FATFS to SPI bridge
-------------------------------------------------------------------------
/* Copyright (c) 2009, Martin Thomas, ChaN
All rights reserved.
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
* Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
* Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in
the documentation and/or other materials provided with the
distribution.
* Neither the name of the copyright holders nor the names of
contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
POSSIBILITY OF SUCH DAMAGE. */
-------------------------------------------------------------------------
FATFS to SDIO bridge
-------------------------------------------------------------------------
/* author: ickle */ - No details given with source code. If you have
anything you'd like to include, pleas e-mail us and let us know.
-------------------------------------------------------------------------
8x8 LCD Font
-------------------------------------------------------------------------
/*
http://forum.osdev.org/viewtopic.php?f=2&t=22033
Created Sunday, May 23, 2010 by Quinn Evans
Renamed and updated Monday 24, 2010
This font (Vincent) is released by me into the public domain. I claim no
copyright, and hereby make this software available to the public for any use,
at any time, free of restrictions, legal or otherwise.
*/
-------------------------------------------------------------------------
CC3000 Host Driver Implementation.
-------------------------------------------------------------------------
Copyright (C) 2011 Texas Instruments Incorporated - http://www.ti.com/
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:
Redistributions of source code must retain the above copyright
notice, this list of conditions and the following disclaimer.
Redistributions in binary form must reproduce the above copyright
notice, this list of conditions and the following disclaimer in the
documentation and/or other materials provided with the
distribution.
Neither the name of Texas Instruments Incorporated nor the names of
its contributors may be used to endorse or promote products derived
from this software without specific prior written permission.
THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
"AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

View File

@@ -0,0 +1,57 @@
_____ _
| __|___ ___ ___ _ _|_|___ ___
| __|_ -| . | _| | | | | . |
|_____|___| _|_| |___|_|_|_|___|
|_|
Copyright 2013 Gordon Williams
http://www.espruino.com
--------------------------------------------------------------
There are a few different binaries in this ZIP file, for different
types of Microcontroller:
espruino_#v##_hystm32_24_ve.bin
- 'HY'STM32F103VET6 ARM with 2.4" LCD display
This is available from eBay
espruino_#v##_hystm32_28_rb.bin
- 'HY'STM32F103RBT6 ARM with 2.8" LCD display
This is available from eBay
espruino_#v##_hystm32_32_vc.bin
- 'HY'STM32F103VCT6 ARM with 3.2" LCD display
This is available from eBay
espruino_#v##_olimexino_stm32.bin
- You will need to overwrite the Maple bootloader to install this.
Espruino is now too large to fit in flash alongside it.
- Olimexino-STM32 Arduino form factor board
- Leaf Labs Maple Arduino form factor board
espruino_#v##_stm32vldiscovery.bin
- STM32VLDISCOVERY board
espruino_#v##_stm32f3discovery.bin
- STM32F3DISCOVERY board
espruino_#v##_stm32f4discovery.bin
- STM32F4DISCOVERY board
espruino_#v##_raspberrypi
- Raspberry Pi executable (just copy it to the device and run it)
NOTE: There is GPIO support (which requires you to run Espruino as root)
however there is no Serial, SPI, OneWire or I2C support at the moment so
you're pretty limited!
For information on devices, and on how to flash these binary files on to
each device, please see our website, http://www.espruino.com
NOTES:
* On the STM32F4DISCOVERY the default USART is USART2 (because
USART1 shares some pins with USB). This means you must connect
serial connections to PA2/PA3 NOT PA9/PA10 as you would for
the STM32VLDISCOVERY.

View File

@@ -0,0 +1 @@
This directory contains auto-generated files

View File

@@ -0,0 +1,7 @@
// auto-generated pin info file
// for board LINUX
#include "jspininfo.h"
const JshPinInfo pinInfo[JSH_PIN_COUNT] = {
};

View File

@@ -0,0 +1,39 @@
// auto-generated pin info file
// for board LINUX
#ifndef __JSPININFO_H_
#define __JSPININFO_H_
#include "jspin.h"
#define JSH_PIN_COUNT 0
#define JSH_PORTA_COUNT 0
#define JSH_PORTB_COUNT 0
#define JSH_PORTC_COUNT 0
#define JSH_PORTD_COUNT 0
#define JSH_PORTE_COUNT 0
#define JSH_PORTF_COUNT 0
#define JSH_PORTG_COUNT 0
#define JSH_PORTH_COUNT 0
#define JSH_PORTA_OFFSET -1
#define JSH_PORTB_OFFSET -1
#define JSH_PORTC_OFFSET -1
#define JSH_PORTD_OFFSET -1
#define JSH_PORTE_OFFSET -1
#define JSH_PORTF_OFFSET -1
#define JSH_PORTG_OFFSET -1
#define JSH_PORTH_OFFSET -1
#define JSH_PININFO_FUNCTIONS 0
typedef struct JshPinInfo {
JsvPinInfoPort port;
JsvPinInfoPin pin;
JsvPinInfoAnalog analog; // TODO: maybe we don't need to store analogs separately
JshPinFunction functions[JSH_PININFO_FUNCTIONS];
} PACKED_FLAGS JshPinInfo;
extern const JshPinInfo pinInfo[JSH_PIN_COUNT];
#endif

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,53 @@
// Automatically generated header file for LINUX
// Generated by scripts/build_platform_config.py
#ifndef _PLATFORM_CONFIG_H
#define _PLATFORM_CONFIG_H
#include <rtthread.h>
#define PC_BOARD_ID "LINUX"
#define PC_BOARD_CHIP "LINUX"
#define PC_BOARD_CHIP_FAMILY "LINUX"
// SYSTICK is the counter that counts up and that we use as the real-time clock
// The smaller this is, the longer we spend in interrupts, but also the more we can sleep!
#define SYSTICK_RANGE 0x1000000 // the Maximum (it is a 24 bit counter) - on Olimexino this is about 0.6 sec
#define SYSTICKS_BEFORE_USB_DISCONNECT 2
#define DEFAULT_BUSY_PIN_INDICATOR (Pin)-1 // no indicator
#define DEFAULT_SLEEP_PIN_INDICATOR (Pin)-1 // no indicator
// When to send the message that the IO buffer is getting full
#define IOBUFFER_XOFF ((TXBUFFERMASK)*6/8)
// When to send the message that we can start receiving again
#define IOBUFFER_XON ((TXBUFFERMASK)*3/8)
#define RAM_TOTAL (-1*1024)
#define FLASH_TOTAL (-1*1024)
#define RESIZABLE_JSVARS // Allocate variables in blocks using malloc
#define USARTS 0
#define SPIS 1
#define I2CS 0
#define ADCS 0
#define DACS 0
#define DEFAULT_CONSOLE_DEVICE EV_USBSERIAL
#define IOBUFFERMASK 31 // (max 255) amount of items in event buffer - events take ~9 bytes each
#define TXBUFFERMASK 31 // (max 255)
// definition to avoid compilation when Pin/platform config is not defined
#define IS_PIN_USED_INTERNALLY(PIN) ((false))
#define IS_PIN_A_LED(PIN) ((false))
#define IS_PIN_A_BUTTON(PIN) ((false))
#endif // _PLATFORM_CONFIG_H

View File

@@ -0,0 +1,233 @@
/*
* This file is part of Espruino, a JavaScript interpreter for Microcontrollers
*
* Copyright (C) 2013 Gordon Williams <gw@pur3.co.uk>
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* ----------------------------------------------------------------------------
* This file is designed to be parsed during the build process
*
* Contains built-in functions for Maths
* ----------------------------------------------------------------------------
*/
#include "jswrap_math.h"
/*JSON{ "type":"class",
"class" : "Math",
"description" : "This is a standard JavaScript class that contains useful Maths routines"
}*/
// -------------------------------------------------------------------- Integer
/*JSON{ "type":"staticmethod",
"class" : "Integer", "name" : "valueOf",
"generate" : "jswrap_integer_valueOf",
"description" : "Given a string containing a single character, return the numeric value of it",
"params" : [ [ "character" ,"JsVar", "A string containing a single character"] ],
"return" : ["int", "The integer value of char"]
}*/
JsVarInt jswrap_integer_valueOf(JsVar *v) {
if (!jsvIsString(v) || jsvGetStringLength(v)!=1)
return 0;
return (int)v->varData.str[0];
}
/*JSON{ "type":"variable", "name" : "NaN",
"generate_full" : "NAN",
"return" : ["float", "Not a Number"]
}*/
// -------------------------------------------------------------------- Double
/*JSON{ "type":"staticmethod",
"class" : "Double", "name" : "doubleToIntBits",
"generate_full" : "*(JsVarInt*)&x",
"description" : " Convert the floating point value given into an integer representing the bits contained in it",
"params" : [ [ "x", "float", "A floating point number"] ],
"return" : ["int", "The integer representation of x"]
}*/
// -------------------------------------------------------------------- Math
/*JSON{ "type":"staticproperty",
"class" : "Math", "name" : "E",
"generate_full" : "2.71828182846",
"return" : ["float", "The value of E - 2.71828182846"]
}*/
/*JSON{ "type":"staticproperty",
"class" : "Math", "name" : "PI",
"generate_full" : "3.14159265359",
"return" : ["float", "The value of PI - 3.14159265359"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "abs",
"generate" : "jswrap_math_abs",
"params" : [ [ "x", "float", "A floating point value"] ],
"return" : ["float", "The absolute value of x (eg, ```Math.abs(2)==2```, but also ```Math.abs(-2)==2```)"]
}*/
JsVarFloat jswrap_math_abs(JsVarFloat x) {
return (x<0)?-x:x;
}
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "acos",
"generate" : "acos",
"params" : [ [ "x", "float", "The value to get the arc cosine of"] ],
"return" : ["float", "The arc cosine of x, between 0 and PI"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "asin",
"generate" : "asin",
"params" : [ [ "x", "float", "The value to get the arc sine of"] ],
"return" : ["float", "The arc sine of x, between -PI/2 and PI/2"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "atan",
"generate" : "atan",
"params" : [ [ "x", "float", "The value to get the arc tangent of"] ],
"return" : ["float", "The arc tangent of x, between -PI/2 and PI/2"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "atan2",
"generate" : "atan2",
"params" : [ [ "y", "float", "The Y-part of the angle to get the arc tangent of"],
[ "x", "float", "The X-part of the angle to get the arc tangent of"] ],
"return" : ["float", "The arctangent of Y/X, between -PI and PI"]
}*/
/* we use sin here, not cos, to try and save a bit of code space */
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "cos",
"generate_full" : "sin(jsvGetFloat(theta) + (3.14159265359/2.0))",
"params" : [ [ "theta", "float", "The angle to get the cosine of"] ],
"return" : ["float", "The cosine of theta"]
}*/
#define DBL_MAX 1.7976931348623157E+308
double fs_fmod(double x, double y)
{
double a, b;
const double c = x;
if (0 > c) {
x = -x;
}
if (0 > y) {
y = -y;
}
if (y != 0 && DBL_MAX >= y && DBL_MAX >= x) {
while (x >= y) {
a = x / 2;
b = y;
while (a >= b) {
b *= 2;
}
x -= b;
}
} else {
x = 0;
}
return 0 > c ? -x : x;
}
double jswrap_math_pow(double x, double y)
{
double p;
if (0 > x && fs_fmod(y, 1) == 0) {
if (fs_fmod(y, 2) == 0) {
p = exp(log(-x) * y);
} else {
p = -exp(log(-x) * y);
}
} else {
if (x != 0 || 0 >= y) {
p = exp(log( x) * y);
} else {
p = 0;
}
}
return p;
}
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "pow",
"generate" : "jswrap_math_pow",
"params" : [ [ "x", "float", "The value to raise to the power"],
[ "y", "float", "The power x should be raised to"] ],
"return" : ["float", "x raised to the power y (x^y)"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "random",
"generate_full" : "(JsVarFloat)rand() / (JsVarFloat)RAND_MAX",
"return" : ["float", "A random number between 0 and 1"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "round",
"generate" : "(JsVarInt)round",
"params" : [ [ "x", "float", "The value to round"] ],
"return" : ["int", "x, rounded to the nearest integer"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "sin",
"generate" : "sin",
"params" : [ [ "theta", "float", "The angle to get the sine of"] ],
"return" : ["float", "The sine of theta"]
}*/
/* we could use the real sqrt - but re-use pow to save on code space */
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "sqrt",
"generate_full" : "jswrap_math_pow(jsvGetFloat(x),0.5)",
"params" : [ [ "x", "float", "The value to take the square root of"] ],
"return" : ["float", "The square root of x"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "ceil",
"generate" : "ceil",
"params" : [ [ "x", "float", "The value to round up"] ],
"return" : ["float", "x, rounded upwards to the nearest integer"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "floor",
"generate" : "floor",
"params" : [ [ "x", "float", "The value to round down"] ],
"return" : ["float", "x, rounded downwards to the nearest integer"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "exp",
"generate" : "exp",
"params" : [ [ "x", "float", "The value raise E to the power of"] ],
"return" : ["float", "E^x"]
}*/
/*JSON{ "type":"staticmethod",
"class" : "Math", "name" : "log",
"generate" : "log",
"params" : [ [ "x", "float", "The value to take the logarithm (base E) root of"] ],
"return" : ["float", "The log (base E) of x"]
}*/
/*JSON{ "type":"staticmethod", "ifndef" : "SAVE_ON_FLASH",
"class" : "Math", "name" : "clip",
"generate" : "jswrap_math_clip",
"description" : "Clip a number to be between min and max (inclusive)",
"params" : [ [ "x", "float", "A floating point value to clip"],
[ "min", "float", "The smallest the value should be"],
[ "max", "float", "The largest the value should be"] ],
"return" : ["float", "The value of x, clipped so as not to be below min or above max."]
}*/
JsVarFloat jswrap_math_clip(JsVarFloat x, JsVarFloat min, JsVarFloat max) {
if (x<min) x=min;
if (x>max) x=max;
return x;
}
/*JSON{ "type":"staticmethod", "ifndef" : "SAVE_ON_FLASH",
"class" : "Math", "name" : "wrap",
"generate" : "wrapAround",
"description" : "Wrap a number around if it is less than 0 or greater than or equal to max. For instance you might do: ```Math.wrap(angleInDegrees, 360)```",
"params" : [ [ "x", "float", "A floating point value to wrap"],
[ "max", "float", "The largest the value should be"] ],
"return" : ["float", "The value of x, wrapped so as not to be below min or above max."]
}*/

View File

@@ -0,0 +1,28 @@
/*
* This file is part of Espruino, a JavaScript interpreter for Microcontrollers
*
* Copyright (C) 2013 Gordon Williams <gw@pur3.co.uk>
*
* This Source Code Form is subject to the terms of the Mozilla Public
* License, v. 2.0. If a copy of the MPL was not distributed with this
* file, You can obtain one at http://mozilla.org/MPL/2.0/.
*
* ----------------------------------------------------------------------------
* Contains built-in functions for Maths
* ----------------------------------------------------------------------------
*/
#include "jsutils.h"
#include "jsvar.h"
#ifdef ARM
#include "mconf.h"
#include "protos.h"
#else
#include <math.h>
#endif
JsVarInt jswrap_integer_valueOf(JsVar *v);
JsVarFloat jswrap_math_abs(JsVarFloat x);
double jswrap_math_pow(double x, double y);
JsVarFloat jswrap_math_clip(JsVarFloat x, JsVarFloat min, JsVarFloat max);

View File

@@ -0,0 +1,32 @@
This suite of C language elementary functions offers support for
not-a-number (NaN) and infinity rules, subnormal numbers, and minus
zero as described by IEEE standard 754 and the Numerical C Extensions
Group (NCEG). For a variety of reasons, many computers cannot take
advantage of these features. You can disable any or all of them by
removing the corresponding preprocessor macros. Check the files
mconf.h and const.c carefully to be sure they are appropriate for your
system.
------------------------------------------
http://www.netlib.org/cephes/readme
Some software in this archive may be from the book _Methods and
Programs for Mathematical Functions_ (Prentice-Hall or Simon & Schuster
International, 1989) or from the Cephes Mathematical Library, a
commercial product. In either event, it is copyrighted by the author.
What you see here may be used freely but it comes with no support or
guarantee.
The two known misprints in the book are repaired here in the
source listings for the gamma function and the incomplete beta
integral.
Stephen L. Moshier
moshier@na-net.ornl.gov

View File

@@ -0,0 +1,167 @@
/* acosh.c
*
* Inverse hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* double x, y, acosh();
*
* y = acosh( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic cosine of argument.
*
* If 1 <= x < 1.5, a rational approximation
*
* sqrt(z) * P(z)/Q(z)
*
* where z = x-1, is used. Otherwise,
*
* acosh(x) = log( x + sqrt( (x-1)(x+1) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC 1,3 30000 4.2e-17 1.1e-17
* IEEE 1,3 30000 4.6e-16 8.7e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* acosh domain |x| < 1 NAN
*
*/
/* acosh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
/* acosh(z) = sqrt(x) * R(x), z = x + 1, interval 0 < x < 0.5 */
#include "mconf.h"
#ifdef UNK
const static double P[] = {
1.18801130533544501356E2,
3.94726656571334401102E3,
3.43989375926195455866E4,
1.08102874834699867335E5,
1.10855947270161294369E5
};
const static double Q[] = {
/* 1.00000000000000000000E0,*/
1.86145380837903397292E2,
4.15352677227719831579E3,
2.97683430363289370382E4,
8.29725251988426222434E4,
7.83869920495893927727E4
};
#endif
#ifdef DEC
static unsigned short P[] = {
0041755,0115055,0144002,0146444,
0043166,0132103,0155150,0150302,
0044006,0057360,0003021,0162753,
0044323,0021557,0175225,0056253,
0044330,0101771,0040046,0006636
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0042072,0022467,0126670,0041232,
0043201,0146066,0152142,0034015,
0043750,0110257,0121165,0026100,
0044242,0007103,0034667,0033173,
0044231,0014576,0175573,0017472
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x59a4,0xb900,0xb345,0x405d,
0x1a18,0x7b4d,0xd688,0x40ae,
0x3cbd,0x00c2,0xcbde,0x40e0,
0xab95,0xff52,0x646d,0x40fa,
0xc1b4,0x2804,0x107f,0x40fb
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x0853,0xf5b7,0x44a6,0x4067,
0x4702,0xda8c,0x3986,0x40b0,
0xa588,0xf44e,0x1215,0x40dd,
0xe6cf,0x6736,0x41c8,0x40f4,
0x63e7,0xdf6f,0x232f,0x40f3
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x405d,0xb345,0xb900,0x59a4,
0x40ae,0xd688,0x7b4d,0x1a18,
0x40e0,0xcbde,0x00c2,0x3cbd,
0x40fa,0x646d,0xff52,0xab95,
0x40fb,0x107f,0x2804,0xc1b4
};
static unsigned short Q[] = {
0x4067,0x44a6,0xf5b7,0x0853,
0x40b0,0x3986,0xda8c,0x4702,
0x40dd,0x1215,0xf44e,0xa588,
0x40f4,0x41c8,0x6736,0xe6cf,
0x40f3,0x232f,0xdf6f,0x63e7,
};
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double log ( double );
extern double sqrt ( double );
#else
double log(), sqrt(), polevl(), p1evl();
#endif
extern double LOGE2, INFINITY, NAN;
double acosh(x)
double x;
{
double a, z;
if( x < 1.0 )
{
mtherr( "acosh", DOMAIN );
return(NAN);
}
if( x > 1.0e8 )
{
#ifdef INFINITIES
if( x == INFINITY )
return( INFINITY );
#endif
return( log(x) + LOGE2 );
}
z = x - 1.0;
if( z < 0.5 )
{
a = sqrt(z) * (polevl(z, P, 4) / p1evl(z, Q, 5) );
return( a );
}
a = sqrt( z*(x+1.0) );
return( log(x + a) );
}

View File

@@ -0,0 +1,324 @@
/* asin.c
*
* Inverse circular sine
*
*
*
* SYNOPSIS:
*
* double x, y, asin();
*
* y = asin( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose sine is x.
*
* A rational function of the form x + x**3 P(x**2)/Q(x**2)
* is used for |x| in the interval [0, 0.5]. If |x| > 0.5 it is
* transformed by the identity
*
* asin(x) = pi/2 - 2 asin( sqrt( (1-x)/2 ) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 40000 2.6e-17 7.1e-18
* IEEE -1, 1 10^6 1.9e-16 5.4e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*
*/
/* acos()
*
* Inverse circular cosine
*
*
*
* SYNOPSIS:
*
* double x, y, acos();
*
* y = acos( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between 0 and pi whose cosine
* is x.
*
* Analytically, acos(x) = pi/2 - asin(x). However if |x| is
* near 1, there is cancellation error in subtracting asin(x)
* from pi/2. Hence if x < -0.5,
*
* acos(x) = pi - 2.0 * asin( sqrt((1+x)/2) );
*
* or if x > +0.5,
*
* acos(x) = 2.0 * asin( sqrt((1-x)/2) ).
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1, 1 50000 3.3e-17 8.2e-18
* IEEE -1, 1 10^6 2.2e-16 6.5e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* asin domain |x| > 1 NAN
*/
/* asin.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
/* arcsin(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.625
Peak relative error = 1.2e-18 */
#if UNK
const static double P[6] = {
4.253011369004428248960E-3,
-6.019598008014123785661E-1,
5.444622390564711410273E0,
-1.626247967210700244449E1,
1.956261983317594739197E1,
-8.198089802484824371615E0,
};
const static double Q[5] = {
/* 1.000000000000000000000E0, */
-1.474091372988853791896E1,
7.049610280856842141659E1,
-1.471791292232726029859E2,
1.395105614657485689735E2,
-4.918853881490881290097E1,
};
#endif
#if DEC
static short P[24] = {
0036213,0056330,0057244,0053234,
0140032,0015011,0114762,0160255,
0040656,0035130,0136121,0067313,
0141202,0014616,0170474,0101731,
0041234,0100076,0151674,0111310,
0141003,0025540,0033165,0077246,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0141153,0155310,0055360,0072530,
0041614,0177001,0027764,0101237,
0142023,0026733,0064653,0133266,
0042013,0101264,0023775,0176351,
0141504,0140420,0050660,0036543,
};
#endif
#if IBMPC
static short P[24] = {
0x8ad3,0x0bd4,0x6b9b,0x3f71,
0x5c16,0x333e,0x4341,0xbfe3,
0x2dd9,0x178a,0xc74b,0x4015,
0x907b,0xde27,0x4331,0xc030,
0x9259,0xda77,0x9007,0x4033,
0xafd5,0x06ce,0x656c,0xc020,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x0eab,0x0b5e,0x7b59,0xc02d,
0x9054,0x25fe,0x9fc0,0x4051,
0x76d7,0x6d35,0x65bb,0xc062,
0xbf9d,0x84ff,0x7056,0x4061,
0x07ac,0x0a36,0x9822,0xc048,
};
#endif
#if MIEEE
static short P[24] = {
0x3f71,0x6b9b,0x0bd4,0x8ad3,
0xbfe3,0x4341,0x333e,0x5c16,
0x4015,0xc74b,0x178a,0x2dd9,
0xc030,0x4331,0xde27,0x907b,
0x4033,0x9007,0xda77,0x9259,
0xc020,0x656c,0x06ce,0xafd5,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc02d,0x7b59,0x0b5e,0x0eab,
0x4051,0x9fc0,0x25fe,0x9054,
0xc062,0x65bb,0x6d35,0x76d7,
0x4061,0x7056,0x84ff,0xbf9d,
0xc048,0x9822,0x0a36,0x07ac,
};
#endif
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x))
0 <= x <= 0.5
Peak relative error = 4.2e-18 */
#if UNK
const static double R[5] = {
2.967721961301243206100E-3,
-5.634242780008963776856E-1,
6.968710824104713396794E0,
-2.556901049652824852289E1,
2.853665548261061424989E1,
};
const static double S[4] = {
/* 1.000000000000000000000E0, */
-2.194779531642920639778E1,
1.470656354026814941758E2,
-3.838770957603691357202E2,
3.424398657913078477438E2,
};
#endif
#if DEC
static short R[20] = {
0036102,0077034,0142164,0174103,
0140020,0036222,0147711,0044173,
0040736,0177655,0153631,0171523,
0141314,0106525,0060015,0055474,
0041344,0045422,0003630,0040344,
};
static short S[16] = {
/* 0040200,0000000,0000000,0000000, */
0141257,0112425,0132772,0166136,
0042023,0010315,0075523,0175020,
0142277,0170104,0126203,0017563,
0042253,0034115,0102662,0022757,
};
#endif
#if IBMPC
static short R[20] = {
0x9f08,0x988e,0x4fc3,0x3f68,
0x290f,0x59f9,0x0792,0xbfe2,
0x3e6a,0xbaf3,0xdff5,0x401b,
0xab68,0xac01,0x91aa,0xc039,
0x081d,0x40f3,0x8962,0x403c,
};
static short S[16] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x5d8c,0xb6bf,0xf2a2,0xc035,
0x7f42,0xaf6a,0x6219,0x4062,
0x63ee,0x9590,0xfe08,0xc077,
0x44be,0xb0b6,0x6709,0x4075,
};
#endif
#if MIEEE
static short R[20] = {
0x3f68,0x4fc3,0x988e,0x9f08,
0xbfe2,0x0792,0x59f9,0x290f,
0x401b,0xdff5,0xbaf3,0x3e6a,
0xc039,0x91aa,0xac01,0xab68,
0x403c,0x8962,0x40f3,0x081d,
};
static short S[16] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0xc035,0xf2a2,0xb6bf,0x5d8c,
0x4062,0x6219,0xaf6a,0x7f42,
0xc077,0xfe08,0x9590,0x63ee,
0x4075,0x6709,0xb0b6,0x44be,
};
#endif
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double sqrt ( double );
double asin ( double );
#else
double sqrt(), polevl(), p1evl();
double asin();
#endif
extern double PIO2, PIO4, NAN;
double asin(x)
double x;
{
double a, p, z, zz;
short sign;
if( x > 0 )
{
sign = 1;
a = x;
}
else
{
sign = -1;
a = -x;
}
if( a > 1.0 )
{
mtherr( "asin", DOMAIN );
return( NAN );
}
if( a > 0.625 )
{
/* arcsin(1-x) = pi/2 - sqrt(2x)(1+R(x)) */
zz = 1.0 - a;
p = zz * polevl( zz, R, 4)/p1evl( zz, S, 4);
zz = sqrt(zz+zz);
z = PIO4 - zz;
zz = zz * p - MOREBITS;
z = z - zz;
z = z + PIO4;
}
else
{
if( a < 1.0e-8 )
{
return(x);
}
zz = a * a;
z = zz * polevl( zz, P, 5)/p1evl( zz, Q, 5);
z = a * z + a;
}
if( sign < 0 )
z = -z;
return(z);
}
double acos(x)
double x;
{
double z;
if( (x < -1.0) || (x > 1.0) )
{
mtherr( "acos", DOMAIN );
return( NAN );
}
if( x > 0.5 )
{
return( 2.0 * asin( sqrt(0.5 - 0.5*x) ) );
}
z = PIO4 - asin(x);
z = z + MOREBITS;
z = z + PIO4;
return( z );
}

View File

@@ -0,0 +1,165 @@
/* asinh.c
*
* Inverse hyperbolic sine
*
*
*
* SYNOPSIS:
*
* double x, y, asinh();
*
* y = asinh( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic sine of argument.
*
* If |x| < 0.5, the function is approximated by a rational
* form x + x**3 P(x)/Q(x). Otherwise,
*
* asinh(x) = log( x + sqrt(1 + x*x) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -3,3 75000 4.6e-17 1.1e-17
* IEEE -1,1 30000 3.7e-16 7.8e-17
* IEEE 1,3 30000 2.5e-16 6.7e-17
*
*/
/* asinh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
const static double P[] = {
-4.33231683752342103572E-3,
-5.91750212056387121207E-1,
-4.37390226194356683570E0,
-9.09030533308377316566E0,
-5.56682227230859640450E0
};
const static double Q[] = {
/* 1.00000000000000000000E0,*/
1.28757002067426453537E1,
4.86042483805291788324E1,
6.95722521337257608734E1,
3.34009336338516356383E1
};
#endif
#ifdef DEC
static unsigned short P[] = {
0136215,0173033,0110410,0105475,
0140027,0076361,0020056,0164520,
0140613,0173401,0160136,0053142,
0141021,0070744,0000503,0176261,
0140662,0021550,0073106,0133351
};
static unsigned short Q[] = {
/* 0040200,0000000,0000000,0000000,*/
0041116,0001336,0034120,0173054,
0041502,0065300,0013144,0021231,
0041613,0022376,0035516,0153063,
0041405,0115216,0054265,0004557
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x1168,0x7221,0xbec3,0xbf71,
0xdd2a,0x2405,0xef9e,0xbfe2,
0xcacc,0x3c0b,0x7ee0,0xc011,
0x7f96,0x8028,0x2e3c,0xc022,
0xd6dd,0x0ec8,0x446d,0xc016
};
static unsigned short Q[] = {
/* 0x0000,0x0000,0x0000,0x3ff0,*/
0x1ec5,0xc70a,0xc05b,0x4029,
0x8453,0x02cc,0x4d58,0x4048,
0xdac6,0xc769,0x649f,0x4051,
0xa12e,0xcb16,0xb351,0x4040
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbf71,0xbec3,0x7221,0x1168,
0xbfe2,0xef9e,0x2405,0xdd2a,
0xc011,0x7ee0,0x3c0b,0xcacc,
0xc022,0x2e3c,0x8028,0x7f96,
0xc016,0x446d,0x0ec8,0xd6dd
};
static unsigned short Q[] = {
0x4029,0xc05b,0xc70a,0x1ec5,
0x4048,0x4d58,0x02cc,0x8453,
0x4051,0x649f,0xc769,0xdac6,
0x4040,0xb351,0xcb16,0xa12e
};
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double sqrt ( double );
extern double log ( double );
#else
double log(), sqrt(), polevl(), p1evl();
#endif
extern double LOGE2, INFINITY;
double asinh(xx)
double xx;
{
double a, z, x;
int sign;
#ifdef MINUSZERO
if( xx == 0.0 )
return(xx);
#endif
if( xx < 0.0 )
{
sign = -1;
x = -xx;
}
else
{
sign = 1;
x = xx;
}
if( x > 1.0e8 )
{
#ifdef INFINITIES
if( x == INFINITY )
return(xx);
#endif
return( sign * (log(x) + LOGE2) );
}
z = x * x;
if( x < 0.5 )
{
a = ( polevl(z, P, 4)/p1evl(z, Q, 4) ) * z;
a = a * x + x;
if( sign < 0 )
a = -a;
return(a);
}
a = sqrt( z + 1.0 );
return( sign * log(x + a) );
}

View File

@@ -0,0 +1,393 @@
/* atan.c
*
* Inverse circular tangent
* (arctangent)
*
*
*
* SYNOPSIS:
*
* double x, y, atan();
*
* y = atan( x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle between -pi/2 and +pi/2 whose tangent
* is x.
*
* Range reduction is from three intervals into the interval
* from zero to 0.66. The approximant uses a rational
* function of degree 4/5 of the form x + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10, 10 50000 2.4e-17 8.3e-18
* IEEE -10, 10 10^6 1.8e-16 5.0e-17
*
*/
/* atan2()
*
* Quadrant correct inverse circular tangent
*
*
*
* SYNOPSIS:
*
* double x, y, z, atan2();
*
* z = atan2( y, x );
*
*
*
* DESCRIPTION:
*
* Returns radian angle whose tangent is y/x.
* Define compile time symbol ANSIC = 1 for ANSI standard,
* range -PI < z <= +PI, args (y,x); else ANSIC = 0 for range
* 0 to 2PI, args (x,y).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -10, 10 10^6 2.5e-16 6.9e-17
* See atan.c.
*
*/
/* atan.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
/* arctan(x) = x + x^3 P(x^2)/Q(x^2)
0 <= x <= 0.66
Peak relative error = 2.6e-18 */
#ifdef UNK
const static double P[5] = {
-8.750608600031904122785E-1,
-1.615753718733365076637E1,
-7.500855792314704667340E1,
-1.228866684490136173410E2,
-6.485021904942025371773E1,
};
const static double Q[5] = {
/* 1.000000000000000000000E0, */
2.485846490142306297962E1,
1.650270098316988542046E2,
4.328810604912902668951E2,
4.853903996359136964868E2,
1.945506571482613964425E2,
};
/* tan( 3*pi/8 ) */
const static double T3P8 = 2.41421356237309504880;
#endif
#ifdef DEC
static short P[20] = {
0140140,0001775,0007671,0026242,
0141201,0041242,0155534,0001715,
0141626,0002141,0132100,0011625,
0141765,0142771,0064055,0150453,
0141601,0131517,0164507,0062164,
};
static short Q[20] = {
/* 0040200,0000000,0000000,0000000, */
0041306,0157042,0154243,0000742,
0042045,0003352,0016707,0150452,
0042330,0070306,0113425,0170730,
0042362,0130770,0116602,0047520,
0042102,0106367,0156753,0013541,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {040432,0101171,0114774,0167462,};
#define T3P8 *(double *)T3P8A
#endif
#ifdef IBMPC
static short P[20] = {
0x2594,0xa1f7,0x007f,0xbfec,
0x807a,0x5b6b,0x2854,0xc030,
0x0273,0x3688,0xc08c,0xc052,
0xba25,0x2d05,0xb8bf,0xc05e,
0xec8e,0xfd28,0x3669,0xc050,
};
static short Q[20] = {
/* 0x0000,0x0000,0x0000,0x3ff0, */
0x603c,0x5b14,0xdbc4,0x4038,
0xfa25,0x43b8,0xa0dd,0x4064,
0xbe3b,0xd2e2,0x0e18,0x407b,
0x49ea,0x13b0,0x563f,0x407e,
0x62ec,0xfbbd,0x519e,0x4068,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {0x9de6,0x333f,0x504f,0x4003};
#define T3P8 *(double *)T3P8A
#endif
#ifdef MIEEE
static short P[20] = {
0xbfec,0x007f,0xa1f7,0x2594,
0xc030,0x2854,0x5b6b,0x807a,
0xc052,0xc08c,0x3688,0x0273,
0xc05e,0xb8bf,0x2d05,0xba25,
0xc050,0x3669,0xfd28,0xec8e,
};
static short Q[20] = {
/* 0x3ff0,0x0000,0x0000,0x0000, */
0x4038,0xdbc4,0x5b14,0x603c,
0x4064,0xa0dd,0x43b8,0xfa25,
0x407b,0x0e18,0xd2e2,0xbe3b,
0x407e,0x563f,0x13b0,0x49ea,
0x4068,0x519e,0xfbbd,0x62ec,
};
/* tan( 3*pi/8 ) = 2.41421356237309504880 */
static unsigned short T3P8A[] = {
0x4003,0x504f,0x333f,0x9de6
};
#define T3P8 *(double *)T3P8A
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double atan ( double );
extern double fabs ( double );
extern int signbit ( double );
extern int isnan ( double );
#else
double polevl(), p1evl(), atan(), fabs();
int signbit(), isnan();
#endif
extern double PI, PIO2, PIO4, INFINITY, NEGZERO, MAXNUM;
/* pi/2 = PIO2 + MOREBITS. */
#ifdef DEC
#define MOREBITS 5.721188726109831840122E-18
#else
#define MOREBITS 6.123233995736765886130E-17
#endif
double atan(x)
double x;
{
double y, z;
short sign, flag;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
#ifdef INFINITIES
if(x == INFINITY)
return(PIO2);
if(x == -INFINITY)
return(-PIO2);
#endif
/* make argument positive and save the sign */
sign = 1;
if( x < 0.0 )
{
sign = -1;
x = -x;
}
/* range reduction */
flag = 0;
if( x > T3P8 )
{
y = PIO2;
flag = 1;
x = -( 1.0/x );
}
else if( x <= 0.66 )
{
y = 0.0;
}
else
{
y = PIO4;
flag = 2;
x = (x-1.0)/(x+1.0);
}
z = x * x;
z = z * polevl( z, P, 4 ) / p1evl( z, Q, 5 );
z = x * z + x;
if( flag == 2 )
z += 0.5 * MOREBITS;
else if( flag == 1 )
z += MOREBITS;
y = y + z;
if( sign < 0 )
y = -y;
return(y);
}
/* atan2 */
#ifdef ANSIC
double atan2( y, x )
#else
double atan2( x, y )
#endif
double x, y;
{
double z, w;
short code;
code = 0;
#ifdef NANS
if( isnan(x) )
return(x);
if( isnan(y) )
return(y);
#endif
#ifdef MINUSZERO
if( y == 0.0 )
{
if( signbit(y) )
{
if( x > 0.0 )
z = y;
else if( x < 0.0 )
z = -PI;
else
{
if( signbit(x) )
z = -PI;
else
z = y;
}
}
else /* y is +0 */
{
if( x == 0.0 )
{
if( signbit(x) )
z = PI;
else
z = 0.0;
}
else if( x > 0.0 )
z = 0.0;
else
z = PI;
}
return z;
}
if( x == 0.0 )
{
if( y > 0.0 )
z = PIO2;
else
z = -PIO2;
return z;
}
#endif /* MINUSZERO */
#ifdef INFINITIES
if( x == INFINITY )
{
if( y == INFINITY )
z = 0.25 * PI;
else if( y == -INFINITY )
z = -0.25 * PI;
else if( y < 0.0 )
z = NEGZERO;
else
z = 0.0;
return z;
}
if( x == -INFINITY )
{
if( y == INFINITY )
z = 0.75 * PI;
else if( y <= -INFINITY )
z = -0.75 * PI;
else if( y >= 0.0 )
z = PI;
else
z = -PI;
return z;
}
if( y == INFINITY )
return( PIO2 );
if( y == -INFINITY )
return( -PIO2 );
#endif
if( x < 0.0 )
code = 2;
if( y < 0.0 )
code |= 1;
#ifdef INFINITIES
if( x == 0.0 )
#else
if( fabs(x) <= (fabs(y) / MAXNUM) )
#endif
{
if( code & 1 )
{
#if ANSIC
return( -PIO2 );
#else
return( 3.0*PIO2 );
#endif
}
if( y == 0.0 )
return( 0.0 );
return( PIO2 );
}
if( y == 0.0 )
{
if( code & 2 )
return( PI );
return( 0.0 );
}
switch( code )
{
#if ANSIC
default:
case 0:
case 1: w = 0.0; break;
case 2: w = PI; break;
case 3: w = -PI; break;
#else
default:
case 0: w = 0.0; break;
case 1: w = 2.0 * PI; break;
case 2:
case 3: w = PI; break;
#endif
}
z = w + atan( y/x );
#ifdef MINUSZERO
if( z == 0.0 && y < 0 )
z = NEGZERO;
#endif
return( z );
}

View File

@@ -0,0 +1,156 @@
/* atanh.c
*
* Inverse hyperbolic tangent
*
*
*
* SYNOPSIS:
*
* double x, y, atanh();
*
* y = atanh( x );
*
*
*
* DESCRIPTION:
*
* Returns inverse hyperbolic tangent of argument in the range
* MINLOG to MAXLOG.
*
* If |x| < 0.5, the rational form x + x**3 P(x)/Q(x) is
* employed. Otherwise,
* atanh(x) = 0.5 * log( (1+x)/(1-x) ).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -1,1 50000 2.4e-17 6.4e-18
* IEEE -1,1 30000 1.9e-16 5.2e-17
*
*/
/* atanh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright (C) 1987, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
const static double P[] = {
-8.54074331929669305196E-1,
1.20426861384072379242E1,
-4.61252884198732692637E1,
6.54566728676544377376E1,
-3.09092539379866942570E1
};
const static double Q[] = {
/* 1.00000000000000000000E0,*/
-1.95638849376911654834E1,
1.08938092147140262656E2,
-2.49839401325893582852E2,
2.52006675691344555838E2,
-9.27277618139601130017E1
};
#endif
#ifdef DEC
static unsigned short P[] = {
0140132,0122235,0105775,0130300,
0041100,0127327,0124407,0034722,
0141470,0100113,0115607,0130535,
0041602,0164721,0003257,0013673,
0141367,0043046,0166673,0045750
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0141234,0101326,0015460,0134564,
0041731,0160115,0116451,0032045,
0142171,0153343,0000532,0167226,
0042174,0000665,0077604,0000310,
0141671,0072235,0031114,0074377
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0xb618,0xb17f,0x5493,0xbfeb,
0xe73a,0xf520,0x15da,0x4028,
0xf62c,0x7370,0x1009,0xc047,
0xe2f7,0x20d5,0x5d3a,0x4050,
0x697d,0xddb7,0xe8c4,0xc03e
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x172f,0xc366,0x905a,0xc033,
0x2685,0xb3a5,0x3c09,0x405b,
0x5dd3,0x602b,0x3adc,0xc06f,
0x8019,0xaff0,0x8036,0x406f,
0x8f20,0xa649,0x2e93,0xc057
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0xbfeb,0x5493,0xb17f,0xb618,
0x4028,0x15da,0xf520,0xe73a,
0xc047,0x1009,0x7370,0xf62c,
0x4050,0x5d3a,0x20d5,0xe2f7,
0xc03e,0xe8c4,0xddb7,0x697d
};
static unsigned short Q[] = {
0xc033,0x905a,0xc366,0x172f,
0x405b,0x3c09,0xb3a5,0x2685,
0xc06f,0x3adc,0x602b,0x5dd3,
0x406f,0x8036,0xaff0,0x8019,
0xc057,0x2e93,0xa649,0x8f20
};
#endif
#ifdef ANSIPROT
extern double fabs ( double );
extern double log ( double x );
extern double polevl ( double x, void *P, int N );
extern double p1evl ( double x, void *P, int N );
#else
double fabs(), log(), polevl(), p1evl();
#endif
extern double INFINITY, NAN;
double atanh(x)
double x;
{
double s, z;
#ifdef MINUSZERO
if( x == 0.0 )
return(x);
#endif
z = fabs(x);
if( z >= 1.0 )
{
if( x == 1.0 )
return( INFINITY );
if( x == -1.0 )
return( -INFINITY );
mtherr( "atanh", DOMAIN );
return( NAN );
}
if( z < 1.0e-7 )
return(x);
if( z < 0.5 )
{
z = x * x;
s = x + x * z * (polevl(z, P, 4) / p1evl(z, Q, 5));
return(s);
}
return( 0.5 * log((1.0+x)/(1.0-x)) );
}

View File

@@ -0,0 +1,142 @@
/* cbrt.c
*
* Cube root
*
*
*
* SYNOPSIS:
*
* double x, y, cbrt();
*
* y = cbrt( x );
*
*
*
* DESCRIPTION:
*
* Returns the cube root of the argument, which may be negative.
*
* Range reduction involves determining the power of 2 of
* the argument. A polynomial of degree 2 applied to the
* mantissa, and multiplication by the cube root of 1, 2, or 4
* approximates the root to within about 0.1%. Then Newton's
* iteration is used three times to converge to an accurate
* result.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10,10 200000 1.8e-17 6.2e-18
* IEEE 0,1e308 30000 1.5e-16 5.0e-17
*
*/
/* cbrt.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1991, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
const static double CBRT2 = 1.2599210498948731647672;
const static double CBRT4 = 1.5874010519681994747517;
const static double CBRT2I = 0.79370052598409973737585;
const static double CBRT4I = 0.62996052494743658238361;
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp();
int isnan(), isfinite();
#endif
double cbrt(x)
double x;
{
int e, rem, sign;
double z;
#ifdef NANS
if( isnan(x) )
return x;
#endif
#ifdef INFINITIES
if( !isfinite(x) )
return x;
#endif
if( x == 0 )
return( x );
if( x > 0 )
sign = 1;
else
{
sign = -1;
x = -x;
}
z = x;
/* extract power of 2, leaving
* mantissa between 0.5 and 1
*/
x = frexp( x, &e );
/* Approximate cube root of number between .5 and 1,
* peak relative error = 9.2e-6
*/
x = (((-1.3466110473359520655053e-1 * x
+ 5.4664601366395524503440e-1) * x
- 9.5438224771509446525043e-1) * x
+ 1.1399983354717293273738e0 ) * x
+ 4.0238979564544752126924e-1;
/* exponent divided by 3 */
if( e >= 0 )
{
rem = e;
e /= 3;
rem -= 3*e;
if( rem == 1 )
x *= CBRT2;
else if( rem == 2 )
x *= CBRT4;
}
/* argument less than 1 */
else
{
e = -e;
rem = e;
e /= 3;
rem -= 3*e;
if( rem == 1 )
x *= CBRT2I;
else if( rem == 2 )
x *= CBRT4I;
e = -e;
}
/* multiply by power of 2 */
x = ldexp( x, e );
/* Newton iteration */
x -= ( x - (z/(x*x)) )*0.33333333333333333333;
#ifdef DEC
x -= ( x - (z/(x*x)) )/3.0;
#else
x -= ( x - (z/(x*x)) )*0.33333333333333333333;
#endif
if( sign < 0 )
x = -x;
return(x);
}

View File

@@ -0,0 +1,82 @@
/* chbevl.c
*
* Evaluate Chebyshev series
*
*
*
* SYNOPSIS:
*
* int N;
* double x, y, coef[N], chebevl();
*
* y = chbevl( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates the series
*
* N-1
* - '
* y = > coef[i] T (x/2)
* - i
* i=0
*
* of Chebyshev polynomials Ti at argument x/2.
*
* Coefficients are stored in reverse order, i.e. the zero
* order term is last in the array. Note N is the number of
* coefficients, not the order.
*
* If coefficients are for the interval a to b, x must
* have been transformed to x -> 2(2x - b - a)/(b-a) before
* entering the routine. This maps x from (a, b) to (-1, 1),
* over which the Chebyshev polynomials are defined.
*
* If the coefficients are for the inverted interval, in
* which (a, b) is mapped to (1/b, 1/a), the transformation
* required is x -> 2(2ab/x - b - a)/(b-a). If b is infinity,
* this becomes x -> 4a/x - 1.
*
*
*
* SPEED:
*
* Taking advantage of the recurrence properties of the
* Chebyshev polynomials, the routine requires one more
* addition per loop than evaluating a nested polynomial of
* the same degree.
*
*/
/* chbevl.c */
/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1985, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
double chbevl( x, array, n )
double x;
double array[];
int n;
{
double b0, b1, b2, *p;
int i;
p = array;
b0 = *p++;
b1 = 0.0;
i = n - 1;
do
{
b2 = b1;
b1 = b0;
b0 = x * b1 - b2 + *p++;
}
while( --i );
return( 0.5*(b0-b2) );
}

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,461 @@
/* cmplx.c
*
* Complex number arithmetic
*
*
*
* SYNOPSIS:
*
* typedef struct {
* double r; real part
* double i; imaginary part
* }cmplx;
*
* cmplx *a, *b, *c;
*
* cadd( a, b, c ); c = b + a
* csub( a, b, c ); c = b - a
* cmul( a, b, c ); c = b * a
* cdiv( a, b, c ); c = b / a
* cneg( c ); c = -c
* cmov( b, c ); c = b
*
*
*
* DESCRIPTION:
*
* Addition:
* c.r = b.r + a.r
* c.i = b.i + a.i
*
* Subtraction:
* c.r = b.r - a.r
* c.i = b.i - a.i
*
* Multiplication:
* c.r = b.r * a.r - b.i * a.i
* c.i = b.r * a.i + b.i * a.r
*
* Division:
* d = a.r * a.r + a.i * a.i
* c.r = (b.r * a.r + b.i * a.i)/d
* c.i = (b.i * a.r - b.r * a.i)/d
* ACCURACY:
*
* In DEC arithmetic, the test (1/z) * z = 1 had peak relative
* error 3.1e-17, rms 1.2e-17. The test (y/z) * (z/y) = 1 had
* peak relative error 8.3e-17, rms 2.1e-17.
*
* Tests in the rectangle {-10,+10}:
* Relative error:
* arithmetic function # trials peak rms
* DEC cadd 10000 1.4e-17 3.4e-18
* IEEE cadd 100000 1.1e-16 2.7e-17
* DEC csub 10000 1.4e-17 4.5e-18
* IEEE csub 100000 1.1e-16 3.4e-17
* DEC cmul 3000 2.3e-17 8.7e-18
* IEEE cmul 100000 2.1e-16 6.9e-17
* DEC cdiv 18000 4.9e-17 1.3e-17
* IEEE cdiv 100000 3.7e-16 1.1e-16
*/
/* cmplx.c
* complex number arithmetic
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double fabs ( double );
extern double cabs ( cmplx * );
extern double sqrt ( double );
extern double atan2 ( double, double );
extern double cos ( double );
extern double sin ( double );
extern double sqrt ( double );
extern double frexp ( double, int * );
extern double ldexp ( double, int );
int isnan ( double );
void cdiv ( cmplx *, cmplx *, cmplx * );
void cadd ( cmplx *, cmplx *, cmplx * );
#else
double fabs(), cabs(), sqrt(), atan2(), cos(), sin();
double sqrt(), frexp(), ldexp();
int isnan();
void cdiv(), cadd();
#endif
extern double MAXNUM, MACHEP, PI, PIO2, INFINITY, NAN;
/*
typedef struct
{
double r;
double i;
}cmplx;
*/
cmplx czero = {0.0, 0.0};
extern cmplx czero;
cmplx cone = {1.0, 0.0};
extern cmplx cone;
/* c = b + a */
void cadd( a, b, c )
register cmplx *a, *b;
cmplx *c;
{
c->r = b->r + a->r;
c->i = b->i + a->i;
}
/* c = b - a */
void csub( a, b, c )
register cmplx *a, *b;
cmplx *c;
{
c->r = b->r - a->r;
c->i = b->i - a->i;
}
/* c = b * a */
void cmul( a, b, c )
register cmplx *a, *b;
cmplx *c;
{
double y;
y = b->r * a->r - b->i * a->i;
c->i = b->r * a->i + b->i * a->r;
c->r = y;
}
/* c = b / a */
void cdiv( a, b, c )
register cmplx *a, *b;
cmplx *c;
{
double y, p, q, w;
y = a->r * a->r + a->i * a->i;
p = b->r * a->r + b->i * a->i;
q = b->i * a->r - b->r * a->i;
if( y < 1.0 )
{
w = MAXNUM * y;
if( (fabs(p) > w) || (fabs(q) > w) || (y == 0.0) )
{
c->r = MAXNUM;
c->i = MAXNUM;
mtherr( "cdiv", OVERFLOW );
return;
}
}
c->r = p/y;
c->i = q/y;
}
/* b = a
Caution, a `short' is assumed to be 16 bits wide. */
void cmov( a, b )
void *a, *b;
{
register short *pa, *pb;
int i;
pa = (short *) a;
pb = (short *) b;
i = 8;
do
*pb++ = *pa++;
while( --i );
}
void cneg( a )
register cmplx *a;
{
a->r = -a->r;
a->i = -a->i;
}
/* cabs()
*
* Complex absolute value
*
*
*
* SYNOPSIS:
*
* double cabs();
* cmplx z;
* double a;
*
* a = cabs( &z );
*
*
*
* DESCRIPTION:
*
*
* If z = x + iy
*
* then
*
* a = sqrt( x**2 + y**2 ).
*
* Overflow and underflow are avoided by testing the magnitudes
* of x and y before squaring. If either is outside half of
* the floating point full scale range, both are rescaled.
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -30,+30 30000 3.2e-17 9.2e-18
* IEEE -10,+10 100000 2.7e-16 6.9e-17
*/
/*
Cephes Math Library Release 2.1: January, 1989
Copyright 1984, 1987, 1989 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
/*
typedef struct
{
double r;
double i;
}cmplx;
*/
#ifdef UNK
#define PREC 27
#define MAXEXP 1024
#define MINEXP -1077
#endif
#ifdef DEC
#define PREC 29
#define MAXEXP 128
#define MINEXP -128
#endif
#ifdef IBMPC
#define PREC 27
#define MAXEXP 1024
#define MINEXP -1077
#endif
#ifdef MIEEE
#define PREC 27
#define MAXEXP 1024
#define MINEXP -1077
#endif
double cabs( z )
register cmplx *z;
{
double x, y, b, re, im;
int ex, ey, e;
#ifdef INFINITIES
/* Note, cabs(INFINITY,NAN) = INFINITY. */
if( z->r == INFINITY || z->i == INFINITY
|| z->r == -INFINITY || z->i == -INFINITY )
return( INFINITY );
#endif
#ifdef NANS
if( isnan(z->r) )
return(z->r);
if( isnan(z->i) )
return(z->i);
#endif
re = fabs( z->r );
im = fabs( z->i );
if( re == 0.0 )
return( im );
if( im == 0.0 )
return( re );
/* Get the exponents of the numbers */
x = frexp( re, &ex );
y = frexp( im, &ey );
/* Check if one number is tiny compared to the other */
e = ex - ey;
if( e > PREC )
return( re );
if( e < -PREC )
return( im );
/* Find approximate exponent e of the geometric mean. */
e = (ex + ey) >> 1;
/* Rescale so mean is about 1 */
x = ldexp( re, -e );
y = ldexp( im, -e );
/* Hypotenuse of the right triangle */
b = sqrt( x * x + y * y );
/* Compute the exponent of the answer. */
y = frexp( b, &ey );
ey = e + ey;
/* Check it for overflow and underflow. */
if( ey > MAXEXP )
{
mtherr( "cabs", OVERFLOW );
return( INFINITY );
}
if( ey < MINEXP )
return(0.0);
/* Undo the scaling */
b = ldexp( b, e );
return( b );
}
/* csqrt()
*
* Complex square root
*
*
*
* SYNOPSIS:
*
* void csqrt();
* cmplx z, w;
*
* csqrt( &z, &w );
*
*
*
* DESCRIPTION:
*
*
* If z = x + iy, r = |z|, then
*
* 1/2
* Im w = [ (r - x)/2 ] ,
*
* Re w = y / 2 Im w.
*
*
* Note that -w is also a square root of z. The root chosen
* is always in the upper half plane.
*
* Because of the potential for cancellation error in r - x,
* the result is sharpened by doing a Heron iteration
* (see sqrt.c) in complex arithmetic.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC -10,+10 25000 3.2e-17 9.6e-18
* IEEE -10,+10 100000 3.2e-16 7.7e-17
*
* 2
* Also tested by csqrt( z ) = z, and tested by arguments
* close to the real axis.
*/
void csqrt( z, w )
cmplx *z, *w;
{
cmplx q, s;
double x, y, r, t;
x = z->r;
y = z->i;
if( y == 0.0 )
{
if( x < 0.0 )
{
w->r = 0.0;
w->i = sqrt(-x);
return;
}
else
{
w->r = sqrt(x);
w->i = 0.0;
return;
}
}
if( x == 0.0 )
{
r = fabs(y);
r = sqrt(0.5*r);
if( y > 0 )
w->r = r;
else
w->r = -r;
w->i = r;
return;
}
/* Approximate sqrt(x^2+y^2) - x = y^2/2x - y^4/24x^3 + ... .
* The relative error in the first term is approximately y^2/12x^2 .
*/
if( (fabs(y) < 2.e-4 * fabs(x))
&& (x > 0) )
{
t = 0.25*y*(y/x);
}
else
{
r = cabs(z);
t = 0.5*(r - x);
}
r = sqrt(t);
q.i = r;
q.r = y/(2.0*r);
/* Heron iteration in complex arithmetic */
cdiv( &q, z, &s );
cadd( &q, &s, w );
w->r *= 0.5;
w->i *= 0.5;
}
double hypot( x, y )
double x, y;
{
cmplx z;
z.r = x;
z.i = y;
return( cabs(&z) );
}

View File

@@ -0,0 +1,252 @@
/* const.c
*
* Globally declared constants
*
*
*
* SYNOPSIS:
*
* extern const double nameofconstant;
*
*
*
*
* DESCRIPTION:
*
* This file contains a number of mathematical constants and
* also some needed size parameters of the computer arithmetic.
* The values are supplied as arrays of hexadecimal integers
* for IEEE arithmetic; arrays of octal constants for DEC
* arithmetic; and in a normal decimal scientific notation for
* other machines. The particular notation used is determined
* by a symbol (DEC, IBMPC, or UNK) defined in the include file
* mconf.h.
*
* The default size parameters are as follows.
*
* For DEC and UNK modes:
* MACHEP = 1.38777878078144567553E-17 2**-56
* MAXLOG = 8.8029691931113054295988E1 log(2**127)
* MINLOG = -8.872283911167299960540E1 log(2**-128)
* MAXNUM = 1.701411834604692317316873e38 2**127
*
* For IEEE arithmetic (IBMPC):
* MACHEP = 1.11022302462515654042E-16 2**-53
* MAXLOG = 7.09782712893383996843E2 log(2**1024)
* MINLOG = -7.08396418532264106224E2 log(2**-1022)
* MAXNUM = 1.7976931348623158E308 2**1024
*
* The global symbols for mathematical constants are
* PI = 3.14159265358979323846 pi
* PIO2 = 1.57079632679489661923 pi/2
* PIO4 = 7.85398163397448309616E-1 pi/4
* SQRT2 = 1.41421356237309504880 sqrt(2)
* SQRTH = 7.07106781186547524401E-1 sqrt(2)/2
* LOG2E = 1.4426950408889634073599 1/log(2)
* SQ2OPI = 7.9788456080286535587989E-1 sqrt( 2/pi )
* LOGE2 = 6.93147180559945309417E-1 log(2)
* LOGSQ2 = 3.46573590279972654709E-1 log(2)/2
* THPIO4 = 2.35619449019234492885 3*pi/4
* TWOOPI = 6.36619772367581343075535E-1 2/pi
*
* These lists are subject to change.
*/
/* const.c */
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
#if 1
const double MACHEP = 1.11022302462515654042E-16; /* 2**-53 */
#else
const double MACHEP = 1.38777878078144567553E-17; /* 2**-56 */
#endif
const double UFLOWTHRESH = 2.22507385850720138309E-308; /* 2**-1022 */
#ifdef DENORMAL
const double MAXLOG = 7.09782712893383996732E2; /* log(MAXNUM) */
/* const double MINLOG = -7.44440071921381262314E2; */ /* log(2**-1074) */
const double MINLOG = -7.451332191019412076235E2; /* log(2**-1075) */
#else
const double MAXLOG = 7.08396418532264106224E2; /* log 2**1022 */
const double MINLOG = -7.08396418532264106224E2; /* log 2**-1022 */
#endif
const double MAXNUM = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
const double PI = 3.14159265358979323846; /* pi */
const double PIO2 = 1.57079632679489661923; /* pi/2 */
const double PIO4 = 7.85398163397448309616E-1; /* pi/4 */
const double SQRT2 = 1.41421356237309504880; /* sqrt(2) */
const double SQRTH = 7.07106781186547524401E-1; /* sqrt(2)/2 */
const double LOG2E = 1.4426950408889634073599; /* 1/log(2) */
const double SQ2OPI = 7.9788456080286535587989E-1; /* sqrt( 2/pi ) */
const double LOGE2 = 6.93147180559945309417E-1; /* log(2) */
const double LOGSQ2 = 3.46573590279972654709E-1; /* log(2)/2 */
const double THPIO4 = 2.35619449019234492885; /* 3*pi/4 */
const double TWOOPI = 6.36619772367581343075535E-1; /* 2/pi */
#ifdef INFINITIES
const double INFINITY = 1.0/0.0; /* 99e999; */
#else
const double INFINITY = 1.79769313486231570815E308; /* 2**1024*(1-MACHEP) */
#endif
#ifdef NANS
const double NAN = 1.0/0.0 - 1.0/0.0;
#else
const double NAN = 0.0;
#endif
#ifdef MINUSZERO
const double NEGZERO = -0.0;
#else
const double NEGZERO = 0.0;
#endif
#endif
#ifdef IBMPC
/* 2**-53 = 1.11022302462515654042E-16 */
const unsigned short MACHEP[4] = {0x0000,0x0000,0x0000,0x3ca0};
const unsigned short UFLOWTHRESH[4] = {0x0000,0x0000,0x0000,0x0010};
#ifdef DENORMAL
/* log(MAXNUM) = 7.09782712893383996732224E2 */
const unsigned short MAXLOG[4] = {0x39ef,0xfefa,0x2e42,0x4086};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/*const unsigned short MINLOG[4] = {0x71c3,0x446d,0x4385,0xc087};*/
const unsigned short MINLOG[4] = {0x3052,0xd52d,0x4910,0xc087};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
const unsigned short MAXLOG[4] = {0xbcd2,0xdd7a,0x232b,0x4086};
/* log(2**-1022) = - 7.08396418532264106224E2 */
const unsigned short MINLOG[4] = {0xbcd2,0xdd7a,0x232b,0xc086};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
const unsigned short MAXNUM[4] = {0xffff,0xffff,0xffff,0x7fef};
const unsigned short PI[4] = {0x2d18,0x5444,0x21fb,0x4009};
const unsigned short PIO2[4] = {0x2d18,0x5444,0x21fb,0x3ff9};
const unsigned short PIO4[4] = {0x2d18,0x5444,0x21fb,0x3fe9};
const unsigned short SQRT2[4] = {0x3bcd,0x667f,0xa09e,0x3ff6};
const unsigned short SQRTH[4] = {0x3bcd,0x667f,0xa09e,0x3fe6};
const unsigned short LOG2E[4] = {0x82fe,0x652b,0x1547,0x3ff7};
const unsigned short SQ2OPI[4] = {0x3651,0x33d4,0x8845,0x3fe9};
const unsigned short LOGE2[4] = {0x39ef,0xfefa,0x2e42,0x3fe6};
const unsigned short LOGSQ2[4] = {0x39ef,0xfefa,0x2e42,0x3fd6};
const unsigned short THPIO4[4] = {0x21d2,0x7f33,0xd97c,0x4002};
const unsigned short TWOOPI[4] = {0xc883,0x6dc9,0x5f30,0x3fe4};
#ifdef INFINITIES
const unsigned short INFINITY[4] = {0x0000,0x0000,0x0000,0x7ff0};
#else
const unsigned short INFINITY[4] = {0xffff,0xffff,0xffff,0x7fef};
#endif
#ifdef NANS
const unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x7ffc};
#else
const unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
const unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x8000};
#else
const unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef MIEEE
/* 2**-53 = 1.11022302462515654042E-16 */
const unsigned short MACHEP[4] = {0x3ca0,0x0000,0x0000,0x0000};
const unsigned short UFLOWTHRESH[4] = {0x0010,0x0000,0x0000,0x0000};
#ifdef DENORMAL
/* log(2**1024) = 7.09782712893383996843E2 */
const unsigned short MAXLOG[4] = {0x4086,0x2e42,0xfefa,0x39ef};
/* log(2**-1074) = - -7.44440071921381262314E2 */
/* const unsigned short MINLOG[4] = {0xc087,0x4385,0x446d,0x71c3}; */
const unsigned short MINLOG[4] = {0xc087,0x4910,0xd52d,0x3052};
#else
/* log(2**1022) = 7.08396418532264106224E2 */
const unsigned short MAXLOG[4] = {0x4086,0x232b,0xdd7a,0xbcd2};
/* log(2**-1022) = - 7.08396418532264106224E2 */
const unsigned short MINLOG[4] = {0xc086,0x232b,0xdd7a,0xbcd2};
#endif
/* 2**1024*(1-MACHEP) = 1.7976931348623158E308 */
const unsigned short MAXNUM[4] = {0x7fef,0xffff,0xffff,0xffff};
const unsigned short PI[4] = {0x4009,0x21fb,0x5444,0x2d18};
const unsigned short PIO2[4] = {0x3ff9,0x21fb,0x5444,0x2d18};
const unsigned short PIO4[4] = {0x3fe9,0x21fb,0x5444,0x2d18};
const unsigned short SQRT2[4] = {0x3ff6,0xa09e,0x667f,0x3bcd};
const unsigned short SQRTH[4] = {0x3fe6,0xa09e,0x667f,0x3bcd};
const unsigned short LOG2E[4] = {0x3ff7,0x1547,0x652b,0x82fe};
const unsigned short SQ2OPI[4] = {0x3fe9,0x8845,0x33d4,0x3651};
const unsigned short LOGE2[4] = {0x3fe6,0x2e42,0xfefa,0x39ef};
const unsigned short LOGSQ2[4] = {0x3fd6,0x2e42,0xfefa,0x39ef};
const unsigned short THPIO4[4] = {0x4002,0xd97c,0x7f33,0x21d2};
const unsigned short TWOOPI[4] = {0x3fe4,0x5f30,0x6dc9,0xc883};
#ifdef INFINITIES
const unsigned short INFINITY[4] = {0x7ff0,0x0000,0x0000,0x0000};
#else
const unsigned short INFINITY[4] = {0x7fef,0xffff,0xffff,0xffff};
#endif
#ifdef NANS
const unsigned short NAN[4] = {0x7ff8,0x0000,0x0000,0x0000};
#else
const unsigned short NAN[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#ifdef MINUSZERO
const unsigned short NEGZERO[4] = {0x8000,0x0000,0x0000,0x0000};
#else
const unsigned short NEGZERO[4] = {0x0000,0x0000,0x0000,0x0000};
#endif
#endif
#ifdef DEC
/* 2**-56 = 1.38777878078144567553E-17 */
const unsigned short MACHEP[4] = {0022200,0000000,0000000,0000000};
const unsigned short UFLOWTHRESH[4] = {0x0080,0x0000,0x0000,0x0000};
/* log 2**127 = 88.029691931113054295988 */
const unsigned short MAXLOG[4] = {041660,007463,0143742,025733,};
/* log 2**-128 = -88.72283911167299960540 */
const unsigned short MINLOG[4] = {0141661,071027,0173721,0147572,};
/* 2**127 = 1.701411834604692317316873e38 */
const unsigned short MAXNUM[4] = {077777,0177777,0177777,0177777,};
const unsigned short PI[4] = {040511,007732,0121041,064302,};
const unsigned short PIO2[4] = {040311,007732,0121041,064302,};
const unsigned short PIO4[4] = {040111,007732,0121041,064302,};
const unsigned short SQRT2[4] = {040265,002363,031771,0157145,};
const unsigned short SQRTH[4] = {040065,002363,031771,0157144,};
const unsigned short LOG2E[4] = {040270,0125073,024534,013761,};
const unsigned short SQ2OPI[4] = {040114,041051,0117241,0131204,};
const unsigned short LOGE2[4] = {040061,071027,0173721,0147572,};
const unsigned short LOGSQ2[4] = {037661,071027,0173721,0147572,};
const unsigned short THPIO4[4] = {040426,0145743,0174631,007222,};
const unsigned short TWOOPI[4] = {040042,0174603,067116,042025,};
/* Approximate infinity by MAXNUM. */
const unsigned short INFINITY[4] = {077777,0177777,0177777,0177777,};
const unsigned short NAN[4] = {0000000,0000000,0000000,0000000};
#ifdef MINUSZERO
const unsigned short NEGZERO[4] = {0000000,0000000,0000000,0100000};
#else
const unsigned short NEGZERO[4] = {0000000,0000000,0000000,0000000};
#endif
#endif
#ifndef UNK
extern const unsigned short MACHEP[];
extern const unsigned short UFLOWTHRESH[];
extern const unsigned short MAXLOG[];
extern const unsigned short UNDLOG[];
extern const unsigned short MINLOG[];
extern const unsigned short MAXNUM[];
extern const unsigned short PI[];
extern const unsigned short PIO2[];
extern const unsigned short PIO4[];
extern const unsigned short SQRT2[];
extern const unsigned short SQRTH[];
extern const unsigned short LOG2E[];
extern const unsigned short SQ2OPI[];
extern const unsigned short LOGE2[];
extern const unsigned short LOGSQ2[];
extern const unsigned short THPIO4[];
extern const unsigned short TWOOPI[];
extern const unsigned short INFINITY[];
extern const unsigned short NAN[];
extern const unsigned short NEGZERO[];
#endif

View File

@@ -0,0 +1,83 @@
/* cosh.c
*
* Hyperbolic cosine
*
*
*
* SYNOPSIS:
*
* double x, y, cosh();
*
* y = cosh( x );
*
*
*
* DESCRIPTION:
*
* Returns hyperbolic cosine of argument in the range MINLOG to
* MAXLOG.
*
* cosh(x) = ( exp(x) + exp(-x) )/2.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 4.0e-17 7.7e-18
* IEEE +-MAXLOG 30000 2.6e-16 5.7e-17
*
*
* ERROR MESSAGES:
*
* message condition value returned
* cosh overflow |x| > MAXLOG MAXNUM
*
*
*/
/* cosh.c */
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1985, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double exp ( double );
extern int isnan ( double );
extern int isfinite ( double );
#else
double exp();
int isnan(), isfinite();
#endif
extern double MAXLOG, INFINITY, LOGE2;
double cosh(x)
double x;
{
double y;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x < 0 )
x = -x;
if( x > (MAXLOG + LOGE2) )
{
mtherr( "cosh", OVERFLOW );
return( INFINITY );
}
if( x >= (MAXLOG - LOGE2) )
{
y = exp(0.5 * x);
y = (0.5 * y) * y;
return(y);
}
y = exp(x);
y = 0.5 * (y + 1.0 / y);
return( y );
}

View File

@@ -0,0 +1,99 @@
CFLAGS= /DEBUG/NOLIST
hfiles= mconf.h-
ofiles= acosh.obj-
asin.obj-
asinh.obj-
atan.obj-
atanh.obj-
cbrt.obj-
chbevl.obj-
const.obj-
cosh.obj-
drand.obj-
exp.obj-
exp10.obj-
fabs.obj-
floor.obj-
log.obj-
log10.obj-
polevl.obj-
pow.obj-
powi.obj-
round.obj-
sin.obj-
sinh.obj-
tan.obj-
tanh.obj-
unity.obj-
sqrt.obj-
floor.obj-
polevl.obj-
mtherr.obj
mtst.exe : $(ofiles)
LINK mtst/option
acosh.obj : acosh.c,$(HFILES)
CC $(CFLAGS) acosh
asin.obj : asin.c,$(HFILES)
CC $(CFLAGS) asin
asinh.obj : asinh.c,$(HFILES)
CC $(CFLAGS) asinh
atan.obj : atan.c,$(HFILES)
CC $(CFLAGS) atan
atan.obj : atan.c,$(HFILES)
CC $(CFLAGS) atan
atanh.obj : atanh.c,$(HFILES)
CC $(CFLAGS) atanh
cbrt.obj : cbrt.c,$(HFILES)
CC $(CFLAGS) cbrt
chbevl.obj : chbevl.c,$(HFILES)
CC $(CFLAGS) chbevl
const.obj : const.c,$(HFILES)
CC $(CFLAGS) const
cosh.obj : cosh.c,$(HFILES)
CC $(CFLAGS) cosh
drand.obj : drand.c,$(HFILES)
CC $(CFLAGS) drand
exp.obj : exp.c,$(HFILES)
CC $(CFLAGS) exp
exp10.obj : exp10.c,$(HFILES)
CC $(CFLAGS) exp10
fabs.obj : fabs.c,$(HFILES)
CC $(CFLAGS) fabs
floor.obj : floor.c,$(HFILES)
CC $(CFLAGS) floor
log.obj : log.c,$(HFILES)
CC $(CFLAGS) log
log10.obj : log10.c,$(HFILES)
CC $(CFLAGS) log10
polevl.obj : polevl.c,$(HFILES)
CC $(CFLAGS) polevl
pow.obj : pow.c,$(HFILES)
CC $(CFLAGS) pow
powi.obj : powi.c,$(HFILES)
CC $(CFLAGS) powi
round.obj : round.c,$(HFILES)
CC $(CFLAGS) round
sin.obj : sin.c,$(HFILES)
CC $(CFLAGS) sin
sinh.obj : sinh.c,$(HFILES)
CC $(CFLAGS) sinh
tan.obj : tan.c,$(HFILES)
CC $(CFLAGS) tan
tanh.obj : tanh.c,$(HFILES)
CC $(CFLAGS) tanh
unity.obj : unity.c,$(HFILES)
CC $(CFLAGS) unity
sqrt.obj : sqrt.c,$(HFILES)
CC $(CFLAGS) sqrt
mtherr.obj : mtherr.c,$(HFILES)
CC $(CFLAGS) mtherr

View File

@@ -0,0 +1,161 @@
/* drand.c
*
* Pseudorandom number generator
*
*
*
* SYNOPSIS:
*
* double y, drand();
*
* drand( &y );
*
*
*
* DESCRIPTION:
*
* Yields a random number 1.0 <= y < 2.0.
*
* The three-generator congruential algorithm by Brian
* Wichmann and David Hill (BYTE magazine, March, 1987,
* pp 127-8) is used. The period, given by them, is
* 6953607871644.
*
* Versions invoked by the different arithmetic compile
* time options DEC, IBMPC, and MIEEE, produce
* approximately the same sequences, differing only in the
* least significant bits of the numbers. The UNK option
* implements the algorithm as recommended in the BYTE
* article. It may be used on all computers. However,
* the low order bits of a double precision number may
* not be adequately random, and may vary due to arithmetic
* implementation details on different computers.
*
* The other compile options generate an additional random
* integer that overwrites the low order bits of the double
* precision number. This reduces the period by a factor of
* two but tends to overcome the problems mentioned.
*
*/
/* Three-generator random number algorithm
* of Brian Wichmann and David Hill
* BYTE magazine, March, 1987 pp 127-8
*
* The period, given by them, is (p-1)(q-1)(r-1)/4 = 6.95e12.
*/
#include "mconf.h"
#ifdef ANSIPROT
static int ranwh ( void );
#else
static int ranwh();
#endif
static int sx = 1;
static int sy = 10000;
static int sz = 3000;
static union {
double d;
unsigned short s[4];
} unkans;
/* This function implements the three
* congruential generators.
*/
static int ranwh()
{
int r, s;
/* sx = sx * 171 mod 30269 */
r = sx/177;
s = sx - 177 * r;
sx = 171 * s - 2 * r;
if( sx < 0 )
sx += 30269;
/* sy = sy * 172 mod 30307 */
r = sy/176;
s = sy - 176 * r;
sy = 172 * s - 35 * r;
if( sy < 0 )
sy += 30307;
/* sz = 170 * sz mod 30323 */
r = sz/178;
s = sz - 178 * r;
sz = 170 * s - 63 * r;
if( sz < 0 )
sz += 30323;
/* The results are in static sx, sy, sz. */
return 0;
}
/* drand.c
*
* Random double precision floating point number between 1 and 2.
*
* C callable:
* drand( &x );
*/
int drand( a )
double *a;
{
unsigned short r;
#ifdef DEC
unsigned short s, t;
#endif
/* This algorithm of Wichmann and Hill computes a floating point
* result:
*/
ranwh();
unkans.d = sx/30269.0 + sy/30307.0 + sz/30323.0;
r = unkans.d;
unkans.d -= r;
unkans.d += 1.0;
/* if UNK option, do nothing further.
* Otherwise, make a random 16 bit integer
* to overwrite the least significant word
* of unkans.
*/
#ifdef UNK
/* do nothing */
#else
ranwh();
r = sx * sy + sz;
#endif
#ifdef DEC
/* To make the numbers as similar as possible
* in all arithmetics, the random integer has
* to be inserted 3 bits higher up in a DEC number.
* An alternative would be put it 3 bits lower down
* in all the other number types.
*/
s = unkans.s[2];
t = s & 07; /* save these bits to put in at the bottom */
s &= 0177770;
s |= (r >> 13) & 07;
unkans.s[2] = s;
t |= r << 3;
unkans.s[3] = t;
#endif
#ifdef IBMPC
unkans.s[0] = r;
#endif
#ifdef MIEEE
unkans.s[3] = r;
#endif
*a = unkans.d;
return 0;
}

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,203 @@
/* exp.c
*
* Exponential function
*
*
*
* SYNOPSIS:
*
* double x, y, exp();
*
* y = exp( x );
*
*
*
* DESCRIPTION:
*
* Returns e (2.71828...) raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
*
* x k f
* e = 2 e.
*
* A Pade' form 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
* of degree 2/3 is used to approximate exp(f) in the basic
* interval [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* DEC +- 88 50000 2.8e-17 7.0e-18
* IEEE +- 708 40000 2.0e-16 5.6e-17
*
*
* Error amplification in the exponential function can be
* a serious matter. The error propagation involves
* exp( X(1+delta) ) = exp(X) ( 1 + X*delta + ... ),
* which shows that a 1 lsb error in representing X produces
* a relative error of X times 1 lsb in the function.
* While the routine gives an accurate result for arguments
* that are exactly represented by a double precision
* computer number, the result contains amplified roundoff
* error for large arguments not exactly represented.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < MINLOG 0.0
* exp overflow x > MAXLOG INFINITY
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
/* Exponential function */
#include "mconf.h"
#ifdef UNK
const static double P[] = {
1.26177193074810590878E-4,
3.02994407707441961300E-2,
9.99999999999999999910E-1,
};
const static double Q[] = {
3.00198505138664455042E-6,
2.52448340349684104192E-3,
2.27265548208155028766E-1,
2.00000000000000000009E0,
};
const static double C1 = 6.93145751953125E-1;
const static double C2 = 1.42860682030941723212E-6;
#endif
#ifdef DEC
static unsigned short P[] = {
0035004,0047156,0127442,0057502,
0036770,0033210,0063121,0061764,
0040200,0000000,0000000,0000000,
};
static unsigned short Q[] = {
0033511,0072665,0160662,0176377,
0036045,0070715,0124105,0132777,
0037550,0134114,0142077,0001637,
0040400,0000000,0000000,0000000,
};
static unsigned short sc1[] = {0040061,0071000,0000000,0000000};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0033277,0137216,0075715,0057117};
#define C2 (*(double *)sc2)
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x4be8,0xd5e4,0x89cd,0x3f20,
0x2c7e,0x0cca,0x06d1,0x3f9f,
0x0000,0x0000,0x0000,0x3ff0,
};
static unsigned short Q[] = {
0x5fa0,0xbc36,0x2eb6,0x3ec9,
0xb6c0,0xb508,0xae39,0x3f64,
0xe074,0x9887,0x1709,0x3fcd,
0x0000,0x0000,0x0000,0x4000,
};
static unsigned short sc1[] = {0x0000,0x0000,0x2e40,0x3fe6};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0xabca,0xcf79,0xf7d1,0x3eb7};
#define C2 (*(double *)sc2)
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f20,0x89cd,0xd5e4,0x4be8,
0x3f9f,0x06d1,0x0cca,0x2c7e,
0x3ff0,0x0000,0x0000,0x0000,
};
static unsigned short Q[] = {
0x3ec9,0x2eb6,0xbc36,0x5fa0,
0x3f64,0xae39,0xb508,0xb6c0,
0x3fcd,0x1709,0x9887,0xe074,
0x4000,0x0000,0x0000,0x0000,
};
static unsigned short sc1[] = {0x3fe6,0x2e40,0x0000,0x0000};
#define C1 (*(double *)sc1)
static unsigned short sc2[] = {0x3eb7,0xf7d1,0xcf79,0xabca};
#define C2 (*(double *)sc2)
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), p1evl(), floor(), ldexp();
int isnan(), isfinite();
#endif
extern double LOGE2, LOG2E, MAXLOG, MINLOG, MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
#endif
double exp(x)
double x;
{
double px, xx;
int n;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x > MAXLOG)
{
#ifdef INFINITIES
return( INFINITY );
#else
mtherr( "exp", OVERFLOW );
return( MAXNUM );
#endif
}
if( x < MINLOG )
{
#ifndef INFINITIES
mtherr( "exp", UNDERFLOW );
#endif
return(0.0);
}
/* Express e**x = e**g 2**n
* = e**g e**( n loge(2) )
* = e**( g + n loge(2) )
*/
px = floor( LOG2E * x + 0.5 ); /* floor() truncates toward -infinity. */
n = px;
x -= px * C1;
x -= px * C2;
/* rational approximation for exponential
* of the fractional part:
* e**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * polevl( xx, P, 2 );
x = px/( polevl( xx, Q, 3 ) - px );
x = 1.0 + 2.0 * x;
/* multiply by power of 2 */
x = ldexp( x, n );
return(x);
}

View File

@@ -0,0 +1,223 @@
/* exp10.c
*
* Base 10 exponential function
* (Common antilogarithm)
*
*
*
* SYNOPSIS:
*
* double x, y, exp10();
*
* y = exp10( x );
*
*
*
* DESCRIPTION:
*
* Returns 10 raised to the x power.
*
* Range reduction is accomplished by expressing the argument
* as 10**x = 2**n 10**f, with |f| < 0.5 log10(2).
* The Pade' form
*
* 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*
* is used to approximate 10**f.
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -307,+307 30000 2.2e-16 5.5e-17
* Test result from an earlier version (2.1):
* DEC -38,+38 70000 3.1e-17 7.0e-18
*
* ERROR MESSAGES:
*
* message condition value returned
* exp10 underflow x < -MAXL10 0.0
* exp10 overflow x > MAXL10 MAXNUM
*
* DEC arithmetic: MAXL10 = 38.230809449325611792.
* IEEE arithmetic: MAXL10 = 308.2547155599167.
*
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1991, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
const static double P[] = {
4.09962519798587023075E-2,
1.17452732554344059015E1,
4.06717289936872725516E2,
2.39423741207388267439E3,
};
const static double Q[] = {
/* 1.00000000000000000000E0,*/
8.50936160849306532625E1,
1.27209271178345121210E3,
2.07960819286001865907E3,
};
/* const static double LOG102 = 3.01029995663981195214e-1; */
const static double LOG210 = 3.32192809488736234787e0;
const static double LG102A = 3.01025390625000000000E-1;
const static double LG102B = 4.60503898119521373889E-6;
/* const static double MAXL10 = 38.230809449325611792; */
const static double MAXL10 = 308.2547155599167;
#endif
#ifdef DEC
static unsigned short P[] = {
0037047,0165657,0114061,0067234,
0041073,0166243,0123052,0144643,
0042313,0055720,0024032,0047443,
0043025,0121714,0070232,0050007,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041652,0027756,0071216,0050075,
0042637,0001367,0077263,0136017,
0043001,0174673,0024157,0133416,
};
/*
static unsigned short L102[] = {0037632,0020232,0102373,0147770};
#define LOG102 *(double *)L102
*/
static unsigned short L210[] = {0040524,0115170,0045715,0015613};
#define LOG210 *(double *)L210
static unsigned short L102A[] = {0037632,0020000,0000000,0000000,};
#define LG102A *(double *)L102A
static unsigned short L102B[] = {0033632,0102373,0147767,0114220,};
#define LG102B *(double *)L102B
static unsigned short MXL[] = {0041430,0166131,0047761,0154130,};
#define MAXL10 ( *(double *)MXL )
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x2dd4,0xf306,0xfd75,0x3fa4,
0x5934,0x74c5,0x7d94,0x4027,
0x49e4,0x0503,0x6b7a,0x4079,
0x4a01,0x8e13,0xb479,0x40a2,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xca08,0xce51,0x45fd,0x4055,
0x7782,0xefd6,0xe05e,0x4093,
0xf6e2,0x650d,0x3f37,0x40a0,
};
/*
static unsigned short L102[] = {0x79ff,0x509f,0x4413,0x3fd3};
#define LOG102 *(double *)L102
*/
static unsigned short L210[] = {0xa371,0x0979,0x934f,0x400a};
#define LOG210 *(double *)L210
static unsigned short L102A[] = {0x0000,0x0000,0x4400,0x3fd3,};
#define LG102A *(double *)L102A
static unsigned short L102B[] = {0xf312,0x79fe,0x509f,0x3ed3,};
#define LG102B *(double *)L102B
const static double MAXL10 = 308.2547155599167;
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3fa4,0xfd75,0xf306,0x2dd4,
0x4027,0x7d94,0x74c5,0x5934,
0x4079,0x6b7a,0x0503,0x49e4,
0x40a2,0xb479,0x8e13,0x4a01,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4055,0x45fd,0xce51,0xca08,
0x4093,0xe05e,0xefd6,0x7782,
0x40a0,0x3f37,0x650d,0xf6e2,
};
/*
static unsigned short L102[] = {0x3fd3,0x4413,0x509f,0x79ff};
#define LOG102 *(double *)L102
*/
static unsigned short L210[] = {0x400a,0x934f,0x0979,0xa371};
#define LOG210 *(double *)L210
static unsigned short L102A[] = {0x3fd3,0x4400,0x0000,0x0000,};
#define LG102A *(double *)L102A
static unsigned short L102B[] = {0x3ed3,0x509f,0x79fe,0xf312,};
#define LG102B *(double *)L102B
const static double MAXL10 = 308.2547155599167;
#endif
#ifdef ANSIPROT
extern double floor ( double );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double floor(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
extern double MAXNUM;
#ifdef INFINITIES
extern double INFINITY;
#endif
double exp10(x)
double x;
{
double px, xx;
short n;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x > MAXL10 )
{
#ifdef INFINITIES
return( INFINITY );
#else
mtherr( "exp10", OVERFLOW );
return( MAXNUM );
#endif
}
if( x < -MAXL10 ) /* Would like to use MINLOG but can't */
{
#ifndef INFINITIES
mtherr( "exp10", UNDERFLOW );
#endif
return(0.0);
}
/* Express 10**x = 10**g 2**n
* = 10**g 10**( n log10(2) )
* = 10**( g + n log10(2) )
*/
px = floor( LOG210 * x + 0.5 );
n = px;
x -= px * LG102A;
x -= px * LG102B;
/* rational approximation for exponential
* of the fractional part:
* 10**x = 1 + 2x P(x**2)/( Q(x**2) - P(x**2) )
*/
xx = x * x;
px = x * polevl( xx, P, 3 );
x = px/( p1evl( xx, Q, 3 ) - px );
x = 1.0 + ldexp( x, 1 );
/* multiply by power of 2 */
x = ldexp( x, n );
return(x);
}

View File

@@ -0,0 +1,183 @@
/* exp2.c
*
* Base 2 exponential function
*
*
*
* SYNOPSIS:
*
* double x, y, exp2();
*
* y = exp2( x );
*
*
*
* DESCRIPTION:
*
* Returns 2 raised to the x power.
*
* Range reduction is accomplished by separating the argument
* into an integer k and fraction f such that
* x k f
* 2 = 2 2.
*
* A Pade' form
*
* 1 + 2x P(x**2) / (Q(x**2) - x P(x**2) )
*
* approximates 2**x in the basic range [-0.5, 0.5].
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE -1022,+1024 30000 1.8e-16 5.4e-17
*
*
* See exp.c for comments on error amplification.
*
*
* ERROR MESSAGES:
*
* message condition value returned
* exp underflow x < -MAXL2 0.0
* exp overflow x > MAXL2 MAXNUM
*
* For DEC arithmetic, MAXL2 = 127.
* For IEEE arithmetic, MAXL2 = 1024.
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
const static double P[] = {
2.30933477057345225087E-2,
2.02020656693165307700E1,
1.51390680115615096133E3,
};
const static double Q[] = {
/* 1.00000000000000000000E0,*/
2.33184211722314911771E2,
4.36821166879210612817E3,
};
#define MAXL2 1024.0
#define MINL2 -1024.0
#endif
#ifdef DEC
static unsigned short P[] = {
0036675,0027102,0122327,0053227,
0041241,0116724,0115412,0157355,
0042675,0036404,0101733,0132226,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0042151,0027450,0077732,0160744,
0043210,0100661,0077550,0056560,
};
#define MAXL2 127.0
#define MINL2 -127.0
#endif
#ifdef IBMPC
static unsigned short P[] = {
0xead3,0x549a,0xa5c8,0x3f97,
0x5bde,0x9361,0x33ba,0x4034,
0x7693,0x907b,0xa7a0,0x4097,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x5c3c,0x0ffb,0x25e5,0x406d,
0x0bae,0x2fed,0x1036,0x40b1,
};
#define MAXL2 1024.0
#define MINL2 -1022.0
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f97,0xa5c8,0x549a,0xead3,
0x4034,0x33ba,0x9361,0x5bde,
0x4097,0xa7a0,0x907b,0x7693,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x406d,0x25e5,0x0ffb,0x5c3c,
0x40b1,0x1036,0x2fed,0x0bae,
};
#define MAXL2 1024.0
#define MINL2 -1022.0
#endif
#ifdef ANSIPROT
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern double floor ( double );
extern double ldexp ( double, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double polevl(), p1evl(), floor(), ldexp();
int isnan(), isfinite();
#endif
#ifdef INFINITIES
extern double INFINITY;
#endif
extern double MAXNUM;
double exp2(x)
double x;
{
double px, xx;
short n;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
if( x > MAXL2)
{
#ifdef INFINITIES
return( INFINITY );
#else
mtherr( "exp2", OVERFLOW );
return( MAXNUM );
#endif
}
if( x < MINL2 )
{
#ifndef INFINITIES
mtherr( "exp2", UNDERFLOW );
#endif
return(0.0);
}
xx = x; /* save x */
/* separate into integer and fractional parts */
px = floor(x+0.5);
n = px;
x = x - px;
/* rational approximation
* exp2(x) = 1 + 2xP(xx)/(Q(xx) - P(xx))
* where xx = x**2
*/
xx = x * x;
px = x * polevl( xx, P, 2 );
x = px / ( p1evl( xx, Q, 2 ) - px );
x = 1.0 + ldexp( x, 1 );
/* scale by power of 2 */
x = ldexp( x, n );
return(x);
}

View File

@@ -0,0 +1,56 @@
/* fabs.c
*
* Absolute value
*
*
*
* SYNOPSIS:
*
* double x, y;
*
* y = fabs( x );
*
*
*
* DESCRIPTION:
*
* Returns the absolute value of the argument.
*
*/
#include "mconf.h"
/* Avoid using UNK if possible. */
#ifdef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif
double fabs(x)
double x;
{
union
{
double d;
short i[4];
} u;
u.d = x;
#ifdef IBMPC
u.i[3] &= 0x7fff;
#endif
#ifdef MIEEE
u.i[0] &= 0x7fff;
#endif
#ifdef DEC
u.i[3] &= 0x7fff;
#endif
#ifdef UNK
if( u.d < 0 )
u.d = -u.d;
#endif
return( u.d );
}

View File

@@ -0,0 +1,453 @@
/* ceil()
* floor()
* frexp()
* ldexp()
* signbit()
* isnan()
* isfinite()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* double ceil(), floor(), frexp(), ldexp();
* int signbit(), isnan(), isfinite();
* double x, y;
* int expnt, n;
*
* y = floor(x);
* y = ceil(x);
* y = frexp( x, &expnt );
* y = ldexp( x, n );
* n = signbit(x);
* n = isnan(x);
* n = isfinite(x);
*
*
*
* DESCRIPTION:
*
* All four routines return a double precision floating point
* result.
*
* floor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* ldexp() multiplies x by 2**n.
*
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif
#ifdef DEC
#define EXPMSK 0x807f
#define MEXP 255
#define NBITS 56
#endif
#ifdef IBMPC
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 53
#endif
#ifdef MIEEE
#define EXPMSK 0x800f
#define MEXP 0x7ff
#define NBITS 53
#endif
extern double MAXNUM, NEGZERO;
#ifdef ANSIPROT
double floor ( double );
int isnan ( double );
int isfinite ( double );
double ldexp ( double, int );
#else
double floor();
int isnan(), isfinite();
double ldexp();
#endif
double ceil(x)
double x;
{
double y;
#ifdef UNK
mtherr( "ceil", DOMAIN );
return(0.0);
#endif
#ifdef NANS
if( isnan(x) )
return( x );
#endif
#ifdef INFINITIES
if(!isfinite(x))
return(x);
#endif
y = floor(x);
if( y < x )
y += 1.0;
#ifdef MINUSZERO
if( y == 0.0 && x < 0.0 )
return( NEGZERO );
#endif
return(y);
}
/* Bit clearing masks: */
static unsigned short bmask[] = {
0xffff,
0xfffe,
0xfffc,
0xfff8,
0xfff0,
0xffe0,
0xffc0,
0xff80,
0xff00,
0xfe00,
0xfc00,
0xf800,
0xf000,
0xe000,
0xc000,
0x8000,
0x0000,
};
double floor(x)
double x;
{
union
{
double y;
unsigned short sh[4];
} u;
unsigned short *p;
int e;
#ifdef UNK
mtherr( "floor", DOMAIN );
return(0.0);
#endif
#ifdef NANS
if( isnan(x) )
return( x );
#endif
#ifdef INFINITIES
if(!isfinite(x))
return(x);
#endif
#ifdef MINUSZERO
if(x == 0.0L)
return(x);
#endif
u.y = x;
/* find the exponent (power of 2) */
#ifdef DEC
p = (unsigned short *)&u.sh[0];
e = (( *p >> 7) & 0377) - 0201;
p += 3;
#endif
#ifdef IBMPC
p = (unsigned short *)&u.sh[3];
e = (( *p >> 4) & 0x7ff) - 0x3ff;
p -= 3;
#endif
#ifdef MIEEE
p = (unsigned short *)&u.sh[0];
e = (( *p >> 4) & 0x7ff) - 0x3ff;
p += 3;
#endif
if( e < 0 )
{
if( u.y < 0.0 )
return( -1.0 );
else
return( 0.0 );
}
e = (NBITS -1) - e;
/* clean out 16 bits at a time */
while( e >= 16 )
{
#ifdef IBMPC
*p++ = 0;
#endif
#ifdef DEC
*p-- = 0;
#endif
#ifdef MIEEE
*p-- = 0;
#endif
e -= 16;
}
/* clear the remaining bits */
if( e > 0 )
*p &= bmask[e];
if( (x < 0) && (u.y != x) )
u.y -= 1.0;
return(u.y);
}
double frexp( x, pw2 )
double x;
int *pw2;
{
union
{
double y;
unsigned short sh[4];
} u;
int i;
#ifdef DENORMAL
int k;
#endif
short *q;
u.y = x;
#ifdef UNK
mtherr( "frexp", DOMAIN );
return(0.0);
#endif
#ifdef IBMPC
q = (short *)&u.sh[3];
#endif
#ifdef DEC
q = (short *)&u.sh[0];
#endif
#ifdef MIEEE
q = (short *)&u.sh[0];
#endif
/* find the exponent (power of 2) */
#ifdef DEC
i = ( *q >> 7) & 0377;
if( i == 0 )
{
*pw2 = 0;
return(0.0);
}
i -= 0200;
*pw2 = i;
*q &= 0x807f; /* strip all exponent bits */
*q |= 040000; /* mantissa between 0.5 and 1 */
return(u.y);
#endif
#ifdef IBMPC
i = ( *q >> 4) & 0x7ff;
if( i != 0 )
goto ieeedon;
#endif
#ifdef MIEEE
i = *q >> 4;
i &= 0x7ff;
if( i != 0 )
goto ieeedon;
#ifdef DENORMAL
#else
*pw2 = 0;
return(0.0);
#endif
#endif
#ifndef DEC
/* Number is denormal or zero */
#ifdef DENORMAL
if( u.y == 0.0 )
{
*pw2 = 0;
return( 0.0 );
}
/* Handle denormal number. */
do
{
u.y *= 2.0;
i -= 1;
k = ( *q >> 4) & 0x7ff;
}
while( k == 0 );
i = i + k;
#endif /* DENORMAL */
ieeedon:
i -= 0x3fe;
*pw2 = i;
*q &= 0x800f;
*q |= 0x3fe0;
return( u.y );
#endif
}
double ldexp( x, pw2 )
double x;
int pw2;
{
union
{
double y;
unsigned short sh[4];
} u;
short *q;
int e;
#ifdef UNK
mtherr( "ldexp", DOMAIN );
return(0.0);
#endif
u.y = x;
#ifdef DEC
q = (short *)&u.sh[0];
e = ( *q >> 7) & 0377;
if( e == 0 )
return(0.0);
#else
#ifdef IBMPC
q = (short *)&u.sh[3];
#endif
#ifdef MIEEE
q = (short *)&u.sh[0];
#endif
while( (e = (*q & 0x7ff0) >> 4) == 0 )
{
if( u.y == 0.0 )
{
return( 0.0 );
}
/* Input is denormal. */
if( pw2 > 0 )
{
u.y *= 2.0;
pw2 -= 1;
}
if( pw2 < 0 )
{
if( pw2 < -53 )
return(0.0);
u.y /= 2.0;
pw2 += 1;
}
if( pw2 == 0 )
return(u.y);
}
#endif /* not DEC */
e += pw2;
/* Handle overflow */
#ifdef DEC
if( e > MEXP )
return( MAXNUM );
#else
if( e >= MEXP )
return( 2.0*MAXNUM );
#endif
/* Handle denormalized results */
if( e < 1 )
{
#ifdef DENORMAL
if( e < -53 )
return(0.0);
*q &= 0x800f;
*q |= 0x10;
/* For denormals, significant bits may be lost even
when dividing by 2. Construct 2^-(1-e) so the result
is obtained with only one multiplication. */
u.y *= ldexp(1.0, e-1);
return(u.y);
#else
return(0.0);
#endif
}
else
{
#ifdef DEC
*q &= 0x807f; /* strip all exponent bits */
*q |= (e & 0xff) << 7;
#else
*q &= 0x800f;
*q |= (e & 0x7ff) << 4;
#endif
return(u.y);
}
}

View File

@@ -0,0 +1,289 @@
# MSDOS Microsoft C makefile for Cephes library
CFLAGS=/c
# For large memory model:
#CFLAGS=/c /AL
# Add /FPa to the CFLAGS if you want to use the fast software FPa arithmetic.
#
# Use the following with /FPa if you do not want to use the 80x87 coprocessor
# or software emulator.
#polevl.obj: polevl.c mconf.h
# cl /c /Ox polevl.c
#
# Use the following instead if you want to use an 80x87 chip or
# software emulator for maximum accuracy computation of the
# polynomial expansions:
polevl.obj: polevl.asm mconf.h
masm polevl.asm/r;
floor.obj: floor.asm
masm floor.asm;
#floor.obj: floor.c mconf.h
# cl $(CFLAGS) floor.c
acosh.obj: acosh.c mconf.h
cl $(CFLAGS) acosh.c
airy.obj: airy.c mconf.h
cl $(CFLAGS) airy.c
asin.obj: asin.c mconf.h
cl $(CFLAGS) asin.c
asinh.obj: asinh.c mconf.h
cl $(CFLAGS) asinh.c
atan.obj: atan.c mconf.h
cl $(CFLAGS) atan.c
atanh.obj: atanh.c mconf.h
cl $(CFLAGS) atanh.c
asinh.obj: asinh.c mconf.h
cl $(CFLAGS) asinh.c
bdtr.obj: bdtr.c mconf.h
cl $(CFLAGS) bdtr.c
beta.obj: beta.c mconf.h
cl $(CFLAGS) beta.c
btdtr.obj: btdtr.c mconf.h
cl $(CFLAGS) btdtr.c
cbrt.obj: cbrt.c mconf.h
cl $(CFLAGS) cbrt.c
chbevl.obj: chbevl.c mconf.h
cl $(CFLAGS) chbevl.c
chdtr.obj: chdtr.c mconf.h
cl $(CFLAGS) chdtr.c
clog.obj: clog.c mconf.h
cl $(CFLAGS) clog.c
cmplx.obj: cmplx.c mconf.h
cl $(CFLAGS) cmplx.c
const.obj: const.c mconf.h
cl $(CFLAGS) const.c
cosh.obj: cosh.c mconf.h
cl $(CFLAGS) cosh.c
dawsn.obj: dawsn.c mconf.h
cl $(CFLAGS) dawsn.c
drand.obj: drand.c mconf.h
cl $(CFLAGS) drand.c
ellie.obj: ellie.c mconf.h
cl $(CFLAGS) ellie.c
ellik.obj: ellik.c mconf.h
cl $(CFLAGS) ellik.c
ellpe.obj: ellpe.c mconf.h
cl $(CFLAGS) ellpe.c
ellpj.obj: ellpj.c mconf.h
cl $(CFLAGS) ellpj.c
ellpk.obj: ellpk.c mconf.h
cl $(CFLAGS) ellpk.c
exp.obj: exp.c mconf.h
cl $(CFLAGS) exp.c
exp10.obj: exp10.c mconf.h
cl $(CFLAGS) exp10.c
exp2.obj: exp2.c mconf.h
cl $(CFLAGS) exp2.c
expn.obj: expn.c mconf.h
cl $(CFLAGS) expn.c
fabs.obj: fabs.c mconf.h
cl $(CFLAGS) fabs.c
fac.obj: fac.c mconf.h
cl $(CFLAGS) fac.c
fdtr.obj: fdtr.c mconf.h
cl $(CFLAGS) fdtr.c
fresnl.obj: fresnl.c mconf.h
cl $(CFLAGS) fresnl.c
gamma.obj: gamma.c mconf.h
cl $(CFLAGS) gamma.c
gdtr.obj: gdtr.c mconf.h
cl $(CFLAGS) gdtr.c
hyp2f1.obj: hyp2f1.c mconf.h
cl $(CFLAGS) hyp2f1.c
hyperg.obj: hyperg.c mconf.h
cl $(CFLAGS) hyperg.c
i0.obj: i0.c mconf.h
cl $(CFLAGS) i0.c
i1.obj: i1.c mconf.h
cl $(CFLAGS) i1.c
igam.obj: igam.c mconf.h
cl $(CFLAGS) igam.c
igami.obj: igami.c mconf.h
cl $(CFLAGS) igami.c
incbet.obj: incbet.c mconf.h
cl $(CFLAGS) incbet.c
incbi.obj: incbi.c mconf.h
cl $(CFLAGS) incbi.c
isnan.obj: isnan.c mconf.h
cl $(CFLAGS) isnan.c
iv.obj: iv.c mconf.h
cl $(CFLAGS) iv.c
j0.obj: j0.c mconf.h
cl $(CFLAGS) j0.c
j1.obj: j1.c mconf.h
cl $(CFLAGS) j1.c
jn.obj: jn.c mconf.h
cl $(CFLAGS) jn.c
jv.obj: jv.c mconf.h
cl $(CFLAGS) jv.c
k0.obj: k0.c mconf.h
cl $(CFLAGS) k0.c
k1.obj: k1.c mconf.h
cl $(CFLAGS) k1.c
kn.obj: kn.c mconf.h
cl $(CFLAGS) kn.c
log.obj: log.c mconf.h
cl $(CFLAGS) log.c
log2.obj: log2.c mconf.h
cl $(CFLAGS) log2.c
log10.obj: log10.c mconf.h
cl $(CFLAGS) log10.c
mtherr.obj: mtherr.c mconf.h
cl $(CFLAGS) mtherr.c
nbdtr.obj: nbdtr.c mconf.h
cl $(CFLAGS) nbdtr.c
ndtr.obj: ndtr.c mconf.h
cl $(CFLAGS) ndtr.c
ndtri.obj: ndtri.c mconf.h
cl $(CFLAGS) ndtri.c
pdtr.obj: pdtr.c mconf.h
cl $(CFLAGS) pdtr.c
pow.obj: pow.c mconf.h
cl $(CFLAGS) pow.c
powi.obj: powi.c mconf.h
cl $(CFLAGS) powi.c
psi.obj: psi.c mconf.h
cl $(CFLAGS) psi.c
rgamma.obj: rgamma.c mconf.h
cl $(CFLAGS) rgamma.c
round.obj: round.c mconf.h
cl $(CFLAGS) round.c
setprec.obj: setprec.87
masm setprec.87;
shichi.obj: shichi.c mconf.h
cl $(CFLAGS) shichi.c
sici.obj: sici.c mconf.h
cl $(CFLAGS) sici.c
sin.obj: sin.c mconf.h
cl $(CFLAGS) sin.c
sindg.obj: sindg.c mconf.h
cl $(CFLAGS) sindg.c
sinh.obj: sinh.c mconf.h
cl $(CFLAGS) sinh.c
spence.obj: spence.c mconf.h
cl $(CFLAGS) spence.c
sqrt.obj: sqrt.87
masm sqrt.87;
#sqrt.obj: sqrt.c
# cl $(CFLAGS) sqrt.c
stdtr.obj: stdtr.c mconf.h
cl $(CFLAGS) stdtr.c
struve.obj: struve.c mconf.h
cl $(CFLAGS) struve.c
tan.obj: tan.c mconf.h
cl $(CFLAGS) tan.c
tandg.obj: tandg.c mconf.h
cl $(CFLAGS) tandg.c
tanh.obj: tanh.c mconf.h
cl $(CFLAGS) tanh.c
yn.obj: yn.c mconf.h
cl $(CFLAGS) yn.c
zeta.obj: zeta.c mconf.h
cl $(CFLAGS) zeta.c
zetac.obj: zetac.c mconf.h
cl $(CFLAGS) zetac.c
polyn.obj: polyn.c mconf.h
cl $(CFLAGS) polyn.c
polmisc.obj: polmisc.c mconf.h
cl $(CFLAGS) polmisc.c
unity.obj: unity.c mconf.h
cl $(CFLAGS) unity.c
fti.lib: acosh.obj airy.obj asin.obj asinh.obj atan.obj atanh.obj bdtr.obj \
beta.obj btdtr.obj cbrt.obj chbevl.obj chdtr.obj clog.obj \
cmplx.obj const.obj cosh.obj dawsn.obj drand.obj ellie.obj ellik.obj \
ellpe.obj ellpj.obj ellpk.obj exp.obj exp10.obj \
exp2.obj expn.obj fabs.obj fac.obj fdtr.obj floor.obj fresnl.obj gamma.obj \
gdtr.obj hyp2f1.obj hyperg.obj i0.c i1.c igam.c igami.obj incbet.obj \
incbi.obj isnan.obj iv.obj j0.obj j1.obj jn.obj jv.obj k0.obj k1.obj \
kn.obj log.obj log2.obj log10.obj mtherr.obj nbdtr.obj ndtr.obj ndtri.obj \
pdtr.obj polevl.obj polmisc.obj polyn.obj pow.obj powi.obj psi.obj \
rgamma.obj round.obj shichi.obj sici.obj sin.obj sindg.obj sinh.obj \
spence.obj sqrt.obj stdtr.obj setprec.obj struve.obj tan.obj \
tandg.obj tanh.obj unity.obj yn.obj zeta.obj zetac.obj \
mconf.h
lib @ftilib.rsp

View File

@@ -0,0 +1,17 @@
fti
y
acosh airy asin asinh atan &
atanh bdtr beta btdtr cbrt chbevl &
chdtr clog cmplx const &
cosh dawsn drand ellie ellik ellpe ellpk &
ellpj exp exp2 exp10 expn fac &
fdtr fresnl gamma gdtr &
hyperg hyp2f1 incbet incbi igam igami isnan &
iv i0 i1 jn jv j0 j1 k0 k1 kn log log2 log10 &
mtherr nbdtr ndtr ndtri pdtr &
polmisc polyn pow powi psi &
rgamma round shichi sici sin sindg &
sinh spence sqrt stdtr struve tan tandg &
tanh unity yn zeta zetac floor fabs polevl
fti.lst
fti

View File

@@ -0,0 +1,237 @@
/* isnan()
* signbit()
* isfinite()
*
* Floating point numeric utilities
*
*
*
* SYNOPSIS:
*
* double ceil(), floor(), frexp(), ldexp();
* int signbit(), isnan(), isfinite();
* double x, y;
* int expnt, n;
*
* y = floor(x);
* y = ceil(x);
* y = frexp( x, &expnt );
* y = ldexp( x, n );
* n = signbit(x);
* n = isnan(x);
* n = isfinite(x);
*
*
*
* DESCRIPTION:
*
* All four routines return a double precision floating point
* result.
*
* floor() returns the largest integer less than or equal to x.
* It truncates toward minus infinity.
*
* ceil() returns the smallest integer greater than or equal
* to x. It truncates toward plus infinity.
*
* frexp() extracts the exponent from x. It returns an integer
* power of two to expnt and the significand between 0.5 and 1
* to y. Thus x = y * 2**expn.
*
* ldexp() multiplies x by 2**n.
*
* signbit(x) returns 1 if the sign bit of x is 1, else 0.
*
* These functions are part of the standard C run time library
* for many but not all C compilers. The ones supplied are
* written in C for either DEC or IEEE arithmetic. They should
* be used only if your compiler library does not already have
* them.
*
* The IEEE versions assume that denormal numbers are implemented
* in the arithmetic. Some modifications will be required if
* the arithmetic has abrupt rather than gradual underflow.
*/
/*
Cephes Math Library Release 2.3: March, 1995
Copyright 1984, 1995 by Stephen L. Moshier
*/
#include "mconf.h"
#ifdef UNK
/* ceil(), floor(), frexp(), ldexp() may need to be rewritten. */
#undef UNK
#if BIGENDIAN
#define MIEEE 1
#else
#define IBMPC 1
#endif
#endif
/* Return 1 if the sign bit of x is 1, else 0. */
int signbit(x)
double x;
{
union
{
double d;
short s[4];
int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
return( u.i[1] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.i[0] < 0 );
#endif
}
else
{
#ifdef IBMPC
return( u.s[3] < 0 );
#endif
#ifdef DEC
return( u.s[3] < 0 );
#endif
#ifdef MIEEE
return( u.s[0] < 0 );
#endif
}
}
/* Return 1 if x is a number that is Not a Number, else return 0. */
int isnan(x)
double x;
{
#ifdef NANS
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( ((u.i[1] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[1] & 0x000fffff) != 0) || (u.i[0] != 0)))
return 1;
#endif
#ifdef DEC
if( (u.s[1] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( ((u.i[0] & 0x7ff00000) == 0x7ff00000)
&& (((u.i[0] & 0x000fffff) != 0) || (u.i[1] != 0)))
return 1;
#endif
return(0);
}
else
{ /* size int not 4 */
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[3] & 0x000f) | u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) == 0)
{
if( (u.s[2] | u.s[1] | u.s[0]) != 0 )
return(1);
}
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) == 0x7ff0)
{
if( ((u.s[0] & 0x000f) | u.s[1] | u.s[2] | u.s[3]) != 0 )
return(1);
}
#endif
return(0);
} /* size int not 4 */
#else
/* No NANS. */
return(0);
#endif
}
/* Return 1 if x is not infinite and is not a NaN. */
int isfinite(x)
double x;
{
#ifdef INFINITIES
union
{
double d;
unsigned short s[4];
unsigned int i[2];
} u;
u.d = x;
if( sizeof(int) == 4 )
{
#ifdef IBMPC
if( (u.i[1] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.i[0] & 0x7ff00000) != 0x7ff00000)
return 1;
#endif
return(0);
}
else
{
#ifdef IBMPC
if( (u.s[3] & 0x7ff0) != 0x7ff0)
return 1;
#endif
#ifdef DEC
if( (u.s[3] & 0x7fff) != 0)
return 1;
#endif
#ifdef MIEEE
if( (u.s[0] & 0x7ff0) != 0x7ff0)
return 1;
#endif
return(0);
}
#else
/* No INFINITY. */
return(1);
#endif
}

View File

@@ -0,0 +1,341 @@
/* log.c
*
* Natural logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log();
*
* y = log( x );
*
*
*
* DESCRIPTION:
*
* Returns the base e (2.718...) logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the logarithm
* of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 150000 1.44e-16 5.06e-17
* IEEE +-MAXNUM 30000 1.20e-16 4.78e-17
* DEC 0, 10 170000 1.8e-17 6.3e-18
*
* In the tests over the interval [+-MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log singularity: x = 0; returns -INFINITY
* log domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
static char fname[] = {"log"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
const static double P[] = {
1.01875663804580931796E-4,
4.97494994976747001425E-1,
4.70579119878881725854E0,
1.44989225341610930846E1,
1.79368678507819816313E1,
7.70838733755885391666E0,
};
const static double Q[] = {
/* 1.00000000000000000000E0, */
1.12873587189167450590E1,
4.52279145837532221105E1,
8.29875266912776603211E1,
7.11544750618563894466E1,
2.31251620126765340583E1,
};
#endif
#ifdef DEC
static unsigned short P[] = {
0037777,0127270,0162547,0057274,
0041001,0054665,0164317,0005341,
0041451,0034104,0031640,0105773,
0041677,0011276,0123617,0160135,
0041701,0126603,0053215,0117250,
0041420,0115777,0135206,0030232,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041220,0144332,0045272,0174241,
0041742,0164566,0035720,0130431,
0042246,0126327,0166065,0116357,
0042372,0033420,0157525,0124560,
0042271,0167002,0066537,0172303,
0041730,0164777,0113711,0044407,
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x1bb0,0x93c3,0xb4c2,0x3f1a,
0x52f2,0x3f56,0xd6f5,0x3fdf,
0x6911,0xed92,0xd2ba,0x4012,
0xeb2e,0xc63e,0xff72,0x402c,
0xc84d,0x924b,0xefd6,0x4031,
0xdcf8,0x7d7e,0xd563,0x401e,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xef8e,0xae97,0x9320,0x4026,
0xc033,0x4e19,0x9d2c,0x4046,
0xbdbd,0xa326,0xbf33,0x4054,
0xae21,0xeb5e,0xc9e2,0x4051,
0x25b2,0x9e1f,0x200a,0x4037,
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f1a,0xb4c2,0x93c3,0x1bb0,
0x3fdf,0xd6f5,0x3f56,0x52f2,
0x4012,0xd2ba,0xed92,0x6911,
0x402c,0xff72,0xc63e,0xeb2e,
0x4031,0xefd6,0x924b,0xc84d,
0x401e,0xd563,0x7d7e,0xdcf8,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4026,0x9320,0xae97,0xef8e,
0x4046,0x9d2c,0x4e19,0xc033,
0x4054,0xbf33,0xa326,0xbdbd,
0x4051,0xc9e2,0xeb5e,0xae21,
0x4037,0x200a,0x9e1f,0x25b2,
};
#endif
/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
const static double R[3] = {
-7.89580278884799154124E-1,
1.63866645699558079767E1,
-6.41409952958715622951E1,
};
const static double S[3] = {
/* 1.00000000000000000000E0,*/
-3.56722798256324312549E1,
3.12093766372244180303E2,
-7.69691943550460008604E2,
};
#endif
#ifdef DEC
static unsigned short R[12] = {
0140112,0020756,0161540,0072035,
0041203,0013743,0114023,0155527,
0141600,0044060,0104421,0050400,
};
static unsigned short S[12] = {
/*0040200,0000000,0000000,0000000,*/
0141416,0130152,0017543,0064122,
0042234,0006000,0104527,0020155,
0142500,0066110,0146631,0174731,
};
#endif
#ifdef IBMPC
static unsigned short R[12] = {
0x0e84,0xdc6c,0x443d,0xbfe9,
0x7b6b,0x7302,0x62fc,0x4030,
0x2a20,0x1122,0x0906,0xc050,
};
static unsigned short S[12] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6d0a,0x43ec,0xd60d,0xc041,
0xe40e,0x112a,0x8180,0x4073,
0x3f3b,0x19b3,0x0d89,0xc088,
};
#endif
#ifdef MIEEE
static unsigned short R[12] = {
0xbfe9,0x443d,0xdc6c,0x0e84,
0x4030,0x62fc,0x7302,0x7b6b,
0xc050,0x0906,0x1122,0x2a20,
};
static unsigned short S[12] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0xc041,0xd60d,0x43ec,0x6d0a,
0x4073,0x8180,0x112a,0xe40e,
0xc088,0x0d89,0x19b3,0x3f3b,
};
#endif
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
#define SQRTH 0.70710678118654752440
extern double INFINITY, NAN;
double log(x)
double x;
{
int e;
#ifdef DEC
short *q;
#endif
double y, z;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
/* rational form */
z = x*x;
z = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
y = e;
z = z - y * 2.121944400546905827679e-4;
z = z + x;
z = z + e * 0.693359375;
goto ldone;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
#if DEC
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 6 ) );
#else
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) );
#endif
if( e )
y = y - e * 2.121944400546905827679e-4;
y = y - ldexp( z, -1 ); /* y - 0.5 * z */
z = x + y;
if( e )
z = z + e * 0.693359375;
ldone:
return( z );
}

View File

@@ -0,0 +1,250 @@
/* log10.c
*
* Common logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log10();
*
* y = log10( x );
*
*
*
* DESCRIPTION:
*
* Returns logarithm to the base 10 of x.
*
* The argument is separated into its exponent and fractional
* parts. The logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 1.5e-16 5.0e-17
* IEEE 0, MAXNUM 30000 1.4e-16 4.8e-17
* DEC 1, MAXNUM 50000 2.5e-17 6.0e-18
*
* In the tests over the interval [1, MAXNUM], the logarithms
* of the random arguments were uniformly distributed over
* [0, MAXLOG].
*
* ERROR MESSAGES:
*
* log10 singularity: x = 0; returns -INFINITY
* log10 domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
static char fname[] = {"log10"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
const static double P[] = {
4.58482948458143443514E-5,
4.98531067254050724270E-1,
6.56312093769992875930E0,
2.97877425097986925891E1,
6.06127134467767258030E1,
5.67349287391754285487E1,
1.98892446572874072159E1
};
const static double Q[] = {
/* 1.00000000000000000000E0, */
1.50314182634250003249E1,
8.27410449222435217021E1,
2.20664384982121929218E2,
3.07254189979530058263E2,
2.14955586696422947765E2,
5.96677339718622216300E1
};
#endif
#ifdef DEC
static unsigned short P[] = {
0034500,0046473,0051374,0135174,
0037777,0037566,0145712,0150321,
0040722,0002426,0031543,0123107,
0041356,0046513,0170752,0004346,
0041562,0071553,0023536,0163343,
0041542,0170221,0024316,0114216,
0041237,0016454,0046611,0104602
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041160,0100260,0067736,0102424,
0041645,0075552,0036563,0147072,
0042134,0125025,0021132,0025320,
0042231,0120211,0046030,0103271,
0042126,0172241,0052151,0120426,
0041556,0125702,0072116,0047103
};
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x974f,0x6a5f,0x09a7,0x3f08,
0x5a1a,0xd979,0xe7ee,0x3fdf,
0x74c9,0xc66c,0x40a2,0x401a,
0x411d,0x7e3d,0xc9a9,0x403d,
0xdcdc,0x64eb,0x4e6d,0x404e,
0xd312,0x2519,0x5e12,0x404c,
0x3130,0x89b1,0xe3a5,0x4033
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xd0a2,0x0dfb,0x1016,0x402e,
0x79c7,0x47ae,0xaf6d,0x4054,
0x455a,0xa44b,0x9542,0x406b,
0x10d7,0x2983,0x3411,0x4073,
0x3423,0x2a8d,0xde94,0x406a,
0xc9c8,0x4e89,0xd578,0x404d
};
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f08,0x09a7,0x6a5f,0x974f,
0x3fdf,0xe7ee,0xd979,0x5a1a,
0x401a,0x40a2,0xc66c,0x74c9,
0x403d,0xc9a9,0x7e3d,0x411d,
0x404e,0x4e6d,0x64eb,0xdcdc,
0x404c,0x5e12,0x2519,0xd312,
0x4033,0xe3a5,0x89b1,0x3130
};
static unsigned short Q[] = {
0x402e,0x1016,0x0dfb,0xd0a2,
0x4054,0xaf6d,0x47ae,0x79c7,
0x406b,0x9542,0xa44b,0x455a,
0x4073,0x3411,0x2983,0x10d7,
0x406a,0xde94,0x2a8d,0x3423,
0x404d,0xd578,0x4e89,0xc9c8
};
#endif
#define SQRTH 0.70710678118654752440
#define L102A 3.0078125E-1
#define L102B 2.48745663981195213739E-4
#define L10EA 4.3359375E-1
#define L10EB 7.00731903251827651129E-4
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
extern double LOGE2, SQRT2, INFINITY, NAN;
double log10(x)
double x;
{
VOLATILE double z;
double y;
#ifdef DEC
short *q;
#endif
int e;
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
/* rational form */
z = x*x;
y = x * ( z * polevl( x, P, 6 ) / p1evl( x, Q, 6 ) );
y = y - ldexp( z, -1 ); /* y - 0.5 * x**2 */
/* multiply log of fraction by log10(e)
* and base 2 exponent by log10(2)
*/
z = (x + y) * L10EB; /* accumulate terms in order of size */
z += y * L10EA;
z += x * L10EA;
z += e * L102B;
z += e * L102A;
return( z );
}

View File

@@ -0,0 +1,348 @@
/* log2.c
*
* Base 2 logarithm
*
*
*
* SYNOPSIS:
*
* double x, y, log2();
*
* y = log2( x );
*
*
*
* DESCRIPTION:
*
* Returns the base 2 logarithm of x.
*
* The argument is separated into its exponent and fractional
* parts. If the exponent is between -1 and +1, the base e
* logarithm of the fraction is approximated by
*
* log(1+x) = x - 0.5 x**2 + x**3 P(x)/Q(x).
*
* Otherwise, setting z = 2(x-1)/x+1),
*
* log(x) = z + z**3 P(z)/Q(z).
*
*
*
* ACCURACY:
*
* Relative error:
* arithmetic domain # trials peak rms
* IEEE 0.5, 2.0 30000 2.0e-16 5.5e-17
* IEEE exp(+-700) 40000 1.3e-16 4.6e-17
*
* In the tests over the interval [exp(+-700)], the logarithms
* of the random arguments were uniformly distributed.
*
* ERROR MESSAGES:
*
* log2 singularity: x = 0; returns -INFINITY
* log2 domain: x < 0; returns NAN
*/
/*
Cephes Math Library Release 2.8: June, 2000
Copyright 1984, 1995, 2000 by Stephen L. Moshier
*/
#include "mconf.h"
static char fname[] = {"log2"};
/* Coefficients for log(1+x) = x - x**2/2 + x**3 P(x)/Q(x)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
const static double P[] = {
1.01875663804580931796E-4,
4.97494994976747001425E-1,
4.70579119878881725854E0,
1.44989225341610930846E1,
1.79368678507819816313E1,
7.70838733755885391666E0,
};
const static double Q[] = {
/* 1.00000000000000000000E0, */
1.12873587189167450590E1,
4.52279145837532221105E1,
8.29875266912776603211E1,
7.11544750618563894466E1,
2.31251620126765340583E1,
};
#define LOG2EA 0.44269504088896340735992
#endif
#ifdef DEC
static unsigned short P[] = {
0037777,0127270,0162547,0057274,
0041001,0054665,0164317,0005341,
0041451,0034104,0031640,0105773,
0041677,0011276,0123617,0160135,
0041701,0126603,0053215,0117250,
0041420,0115777,0135206,0030232,
};
static unsigned short Q[] = {
/*0040200,0000000,0000000,0000000,*/
0041220,0144332,0045272,0174241,
0041742,0164566,0035720,0130431,
0042246,0126327,0166065,0116357,
0042372,0033420,0157525,0124560,
0042271,0167002,0066537,0172303,
0041730,0164777,0113711,0044407,
};
static unsigned short L[5] = {0037742,0124354,0122560,0057703};
#define LOG2EA (*(double *)(&L[0]))
#endif
#ifdef IBMPC
static unsigned short P[] = {
0x1bb0,0x93c3,0xb4c2,0x3f1a,
0x52f2,0x3f56,0xd6f5,0x3fdf,
0x6911,0xed92,0xd2ba,0x4012,
0xeb2e,0xc63e,0xff72,0x402c,
0xc84d,0x924b,0xefd6,0x4031,
0xdcf8,0x7d7e,0xd563,0x401e,
};
static unsigned short Q[] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0xef8e,0xae97,0x9320,0x4026,
0xc033,0x4e19,0x9d2c,0x4046,
0xbdbd,0xa326,0xbf33,0x4054,
0xae21,0xeb5e,0xc9e2,0x4051,
0x25b2,0x9e1f,0x200a,0x4037,
};
static unsigned short L[5] = {0x0bf8,0x94ae,0x551d,0x3fdc};
#define LOG2EA (*(double *)(&L[0]))
#endif
#ifdef MIEEE
static unsigned short P[] = {
0x3f1a,0xb4c2,0x93c3,0x1bb0,
0x3fdf,0xd6f5,0x3f56,0x52f2,
0x4012,0xd2ba,0xed92,0x6911,
0x402c,0xff72,0xc63e,0xeb2e,
0x4031,0xefd6,0x924b,0xc84d,
0x401e,0xd563,0x7d7e,0xdcf8,
};
static unsigned short Q[] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0x4026,0x9320,0xae97,0xef8e,
0x4046,0x9d2c,0x4e19,0xc033,
0x4054,0xbf33,0xa326,0xbdbd,
0x4051,0xc9e2,0xeb5e,0xae21,
0x4037,0x200a,0x9e1f,0x25b2,
};
static unsigned short L[5] = {0x3fdc,0x551d,0x94ae,0x0bf8};
#define LOG2EA (*(double *)(&L[0]))
#endif
/* Coefficients for log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/(x+1)
* 1/sqrt(2) <= x < sqrt(2)
*/
#ifdef UNK
const static double R[3] = {
-7.89580278884799154124E-1,
1.63866645699558079767E1,
-6.41409952958715622951E1,
};
const static double S[3] = {
/* 1.00000000000000000000E0,*/
-3.56722798256324312549E1,
3.12093766372244180303E2,
-7.69691943550460008604E2,
};
/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340735992
#endif
#ifdef DEC
static unsigned short R[12] = {
0140112,0020756,0161540,0072035,
0041203,0013743,0114023,0155527,
0141600,0044060,0104421,0050400,
};
static unsigned short S[12] = {
/*0040200,0000000,0000000,0000000,*/
0141416,0130152,0017543,0064122,
0042234,0006000,0104527,0020155,
0142500,0066110,0146631,0174731,
};
/* log2(e) - 1 */
#define LOG2EA 0.44269504088896340735992L
#endif
#ifdef IBMPC
static unsigned short R[12] = {
0x0e84,0xdc6c,0x443d,0xbfe9,
0x7b6b,0x7302,0x62fc,0x4030,
0x2a20,0x1122,0x0906,0xc050,
};
static unsigned short S[12] = {
/*0x0000,0x0000,0x0000,0x3ff0,*/
0x6d0a,0x43ec,0xd60d,0xc041,
0xe40e,0x112a,0x8180,0x4073,
0x3f3b,0x19b3,0x0d89,0xc088,
};
#endif
#ifdef MIEEE
static unsigned short R[12] = {
0xbfe9,0x443d,0xdc6c,0x0e84,
0x4030,0x62fc,0x7302,0x7b6b,
0xc050,0x0906,0x1122,0x2a20,
};
static unsigned short S[12] = {
/*0x3ff0,0x0000,0x0000,0x0000,*/
0xc041,0xd60d,0x43ec,0x6d0a,
0x4073,0x8180,0x112a,0xe40e,
0xc088,0x0d89,0x19b3,0x3f3b,
};
#endif
#ifdef ANSIPROT
extern double frexp ( double, int * );
extern double ldexp ( double, int );
extern double polevl ( double, void *, int );
extern double p1evl ( double, void *, int );
extern int isnan ( double );
extern int isfinite ( double );
#else
double frexp(), ldexp(), polevl(), p1evl();
int isnan(), isfinite();
#endif
#define SQRTH 0.70710678118654752440
extern double LOGE2, INFINITY, NAN;
double log2(x)
double x;
{
int e;
double y;
VOLATILE double z;
#ifdef DEC
short *q;
#endif
#ifdef NANS
if( isnan(x) )
return(x);
#endif
#ifdef INFINITIES
if( x == INFINITY )
return(x);
#endif
/* Test for domain */
if( x <= 0.0 )
{
if( x == 0.0 )
{
mtherr( fname, SING );
return( -INFINITY );
}
else
{
mtherr( fname, DOMAIN );
return( NAN );
}
}
/* separate mantissa from exponent */
#ifdef DEC
q = (short *)&x;
e = *q; /* short containing exponent */
e = ((e >> 7) & 0377) - 0200; /* the exponent */
*q &= 0177; /* strip exponent from x */
*q |= 040000; /* x now between 0.5 and 1 */
#endif
/* Note, frexp is used so that denormal numbers
* will be handled properly.
*/
#ifdef IBMPC
x = frexp( x, &e );
/*
q = (short *)&x;
q += 3;
e = *q;
e = ((e >> 4) & 0x0fff) - 0x3fe;
*q &= 0x0f;
*q |= 0x3fe0;
*/
#endif
/* Equivalent C language standard library function: */
#ifdef UNK
x = frexp( x, &e );
#endif
#ifdef MIEEE
x = frexp( x, &e );
#endif
/* logarithm using log(x) = z + z**3 P(z)/Q(z),
* where z = 2(x-1)/x+1)
*/
if( (e > 2) || (e < -2) )
{
if( x < SQRTH )
{ /* 2( 2x-1 )/( 2x+1 ) */
e -= 1;
z = x - 0.5;
y = 0.5 * z + 0.5;
}
else
{ /* 2 (x-1)/(x+1) */
z = x - 0.5;
z -= 0.5;
y = 0.5 * x + 0.5;
}
x = z / y;
z = x*x;
y = x * ( z * polevl( z, R, 2 ) / p1evl( z, S, 3 ) );
goto ldone;
}
/* logarithm using log(1+x) = x - .5x**2 + x**3 P(x)/Q(x) */
if( x < SQRTH )
{
e -= 1;
x = ldexp( x, 1 ) - 1.0; /* 2x - 1 */
}
else
{
x = x - 1.0;
}
z = x*x;
#if DEC
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 6 ) ) - ldexp( z, -1 );
#else
y = x * ( z * polevl( x, P, 5 ) / p1evl( x, Q, 5 ) ) - ldexp( z, -1 );
#endif
ldone:
/* Multiply log of fraction by log2(e)
* and base 2 exponent by 1
*
* ***CAUTION***
*
* This sequence of operations is critical and it may
* be horribly defeated by some compiler optimizers.
*/
z = y * LOG2EA;
z += x * LOG2EA;
z += y;
z += x;
z += e;
return( z );
}

View File

@@ -0,0 +1,199 @@
/* mconf.h
*
* Common include file for math routines
*
*
*
* SYNOPSIS:
*
* #include "mconf.h"
*
*
*
* DESCRIPTION:
*
* This file contains definitions for error codes that are
* passed to the common error handling routine mtherr()
* (which see).
*
* The file also includes a conditional assembly definition
* for the type of computer arithmetic (IEEE, DEC, Motorola
* IEEE, or UNKnown).
*
* For Digital Equipment PDP-11 and VAX computers, certain
* IBM systems, and others that use numbers with a 56-bit
* significand, the symbol DEC should be defined. In this
* mode, most floating point constants are given as arrays
* of octal integers to eliminate decimal to binary conversion
* errors that might be introduced by the compiler.
*
* For little-endian computers, such as IBM PC, that follow the
* IEEE Standard for Binary Floating Point Arithmetic (ANSI/IEEE
* Std 754-1985), the symbol IBMPC should be defined. These
* numbers have 53-bit significands. In this mode, constants
* are provided as arrays of hexadecimal 16 bit integers.
*
* Big-endian IEEE format is denoted MIEEE. On some RISC
* systems such as Sun SPARC, double precision constants
* must be stored on 8-byte address boundaries. Since integer
* arrays may be aligned differently, the MIEEE configuration
* may fail on such machines.
*
* To accommodate other types of computer arithmetic, all
* constants are also provided in a normal decimal radix
* which one can hope are correctly converted to a suitable
* format by the available C language compiler. To invoke
* this mode, define the symbol UNK.
*
* An important difference among these modes is a predefined
* set of machine arithmetic constants for each. The numbers
* MACHEP (the machine roundoff error), MAXNUM (largest number
* represented), and several other parameters are preset by
* the configuration symbol. Check the file const.c to
* ensure that these values are correct for your computer.
*
* Configurations NANS, INFINITIES, MINUSZERO, and DENORMAL
* may fail on many systems. Verify that they are supposed
* to work on your computer.
*/
/*
Cephes Math Library Release 2.3: June, 1995
Copyright 1984, 1987, 1989, 1995 by Stephen L. Moshier
*/
/* Define if the `long double' type works. */
#define HAVE_LONG_DOUBLE 1
/* Define as the return type of signal handlers (int or void). */
#define RETSIGTYPE void
/* Define if you have the ANSI C header files. */
#define STDC_HEADERS 1
/* Define if your processor stores words with the most significant
byte first (like Motorola and SPARC, unlike Intel and VAX). */
/* #undef WORDS_BIGENDIAN */
/* Define if floating point words are bigendian. */
/* #undef FLOAT_WORDS_BIGENDIAN */
/* The number of bytes in a int. */
#define SIZEOF_INT 4
/* Define if you have the <string.h> header file. */
#define HAVE_STRING_H 1
/* Name of package */
#define PACKAGE "cephes"
/* Version number of package */
#define VERSION "2.7"
/* Constant definitions for math error conditions
*/
#define DOMAIN 1 /* argument domain error */
#define SING 2 /* argument singularity */
#define OVERFLOW 3 /* overflow range error */
#define UNDERFLOW 4 /* underflow range error */
#define TLOSS 5 /* total loss of precision */
#define PLOSS 6 /* partial loss of precision */
#define EDOM 33
#define ERANGE 34
/* Complex numeral. */
typedef struct
{
double r;
double i;
} cmplx;
#ifdef HAVE_LONG_DOUBLE
/* Long double complex numeral. */
typedef struct
{
long double r;
long double i;
} cmplxl;
#endif
/* Type of computer arithmetic */
/* PDP-11, Pro350, VAX:
*/
/* #define DEC 1 */
/* Intel IEEE, low order words come first:
*/
/* #define IBMPC 1 */
/* Motorola IEEE, high order words come first
* (Sun 680x0 workstation):
*/
/* #define MIEEE 1 */
/* UNKnown arithmetic, invokes coefficients given in
* normal decimal format. Beware of range boundary
* problems (MACHEP, MAXLOG, etc. in const.c) and
* roundoff problems in pow.c:
* (Sun SPARCstation)
*/
#define UNK 1
/* If you define UNK, then be sure to set BIGENDIAN properly. */
#ifdef FLOAT_WORDS_BIGENDIAN
#define BIGENDIAN 1
#else
#define BIGENDIAN 0
#endif
/* Define this `volatile' if your compiler thinks
* that floating point arithmetic obeys the associative
* and distributive laws. It will defeat some optimizations
* (but probably not enough of them).
*
* #define VOLATILE volatile
*/
#define VOLATILE
/* For 12-byte long doubles on an i386, pad a 16-bit short 0
* to the end of real constants initialized by integer arrays.
*
* #define XPD 0,
*
* Otherwise, the type is 10 bytes long and XPD should be
* defined blank (e.g., Microsoft C).
*
* #define XPD
*/
#define XPD 0,
/* Define to support tiny denormal numbers, else undefine. */
#define DENORMAL 1
/* Define to ask for infinity support, else undefine. */
#define INFINITIES 1
/* Define to ask for support of numbers that are Not-a-Number,
else undefine. This may automatically define INFINITIES in some files. */
#define NANS 1
/* Define to distinguish between -0.0 and +0.0. */
#define MINUSZERO 1
/* Define 1 for ANSI C atan2() function
See atan.c and clog.c. */
#define ANSIC 1
/* Get ANSI function prototypes, if you want them. */
#if 1
/* #ifdef __STDC__ */
#define ANSIPROT 1
int mtherr ( char *, int );
#else
int mtherr();
#endif
/* Variable for error reporting. See mtherr.c. */
extern int merror;

View File

@@ -0,0 +1,122 @@
/* Program to test range reduction of trigonometry functions
*
* -- Steve Moshier
*/
#include "mconf.h"
#ifdef ANSIPROT
extern double floor ( double );
extern double ldexp ( double, int );
extern double sin ( double );
#else
double floor(), ldexp(), sin();
#endif
#define TPI 6.283185307179586476925
main()
{
char s[40];
double a, n, t, x, y, z;
int lflg;
x = TPI/4.0;
t = 1.0;
loop:
t = 2.0 * t;
/* Stop testing at a point beyond which the integer part of
* x/2pi cannot be represented exactly by a double precision number.
* The library trigonometry functions will probably give up long before
* this point is reached.
*/
if( t > 1.0e16 )
exit(0);
/* Adjust the following to choose a nontrivial x
* where test function(x) has a slope of about 1 or more.
*/
x = TPI * t + 0.5;
z = x;
lflg = 0;
inlup:
/* floor() returns the largest integer less than its argument.
* If you do not have this, or AINT(), then you may convert x/TPI
* to a long integer and then back to double; but in that case
* x will be limited to the largest value that will fit into a
* long integer.
*/
n = floor( z/TPI );
/* Carefully subtract 2 pi n from x.
* This is done by subtracting n * 2**k in such a way that there
* is no arithmetic cancellation error at any step. The k are the
* bits in the number 2 pi.
*
* If you do not have ldexp(), then you may multiply or
* divide n by an appropriate power of 2 after each step.
* For example:
* a = z - 4*n;
* a -= 2*n;
* n /= 4;
* a -= n; n/4
* n /= 8;
* a -= n; n/32
* etc.
* This will only work if division by a power of 2 is exact.
*/
a = z - ldexp(n, 2); /* 4n */
a -= ldexp( n, 1); /* 2n */
a -= ldexp( n, -2 ); /* n/4 */
a -= ldexp( n, -5 ); /* n/32 */
a -= ldexp( n, -9 ); /* n/512 */
a += ldexp( n, -15 ); /* add n/32768 */
a -= ldexp( n, -17 ); /* n/131072 */
a -= ldexp( n, -18 );
a -= ldexp( n, -20 );
a -= ldexp( n, -22 );
a -= ldexp( n, -24 );
a -= ldexp( n, -28 );
a -= ldexp( n, -32 );
a -= ldexp( n, -37 );
a -= ldexp( n, -39 );
a -= ldexp( n, -40 );
a -= ldexp( n, -42 );
a -= ldexp( n, -46 );
a -= ldexp( n, -47 );
/* Subtract what is left of 2 pi n after all the above reductions.
*/
a -= 2.44929359829470635445e-16 * n;
/* If the test is extended too far, it is possible
* to have chosen the wrong value of n. The following
* will fix that, but at some reduction in accuracy.
*/
if( (a > TPI) || (a < -1e-11) )
{
z = a;
lflg += 1;
printf( "Warning! Reduction failed on first try.\n" );
goto inlup;
}
if( a < 0.0 )
{
printf( "Warning! Reduced value < 0\n" );
a += TPI;
}
/* Compute the test function at x and at a = x mod 2 pi.
*/
y = sin(x);
z = sin(a);
printf( "sin(%.15e) error = %.3e\n", x, y-z );
goto loop;
}

View File

@@ -0,0 +1,103 @@
/* mtherr.c
*
* Library common error handling routine
*
*
*
* SYNOPSIS:
*
* char *fctnam;
* int code;
* int mtherr();
*
* mtherr( fctnam, code );
*
*
*
* DESCRIPTION:
*
* This routine may be called to report one of the following
* error conditions (in the include file mconf.h).
*
* Mnemonic Value Significance
*
* DOMAIN 1 argument domain error
* SING 2 function singularity
* OVERFLOW 3 overflow range error
* UNDERFLOW 4 underflow range error
* TLOSS 5 total loss of precision
* PLOSS 6 partial loss of precision
* EDOM 33 Unix domain error code
* ERANGE 34 Unix range error code
*
* The default version of the file prints the function name,
* passed to it by the pointer fctnam, followed by the
* error condition. The display is directed to the standard
* output device. The routine then returns to the calling
* program. Users may wish to modify the program to abort by
* calling exit() under severe error conditions such as domain
* errors.
*
* Since all error conditions pass control to this function,
* the display may be easily changed, eliminated, or directed
* to an error logging device.
*
* SEE ALSO:
*
* mconf.h
*
*/
/*
Cephes Math Library Release 2.0: April, 1987
Copyright 1984, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
#if 0
#include <stdio.h>
#endif
#include "mconf.h"
int merror = 0;
/* Notice: the order of appearance of the following
* messages is bound to the error codes defined
* in mconf.h.
*/
static char *ermsg[7] = {
"unknown", /* error code 0 */
"domain", /* error code 1 */
"singularity", /* et seq. */
"overflow",
"underflow",
"total loss of precision",
"partial loss of precision"
};
int mtherr( name, code )
char *name;
int code;
{
#if 0
/* Display string passed by calling program,
* which is supposed to be the name of the
* function in which the error occurred:
*/
printf( "\n%s ", name );
/* Set global error message word */
merror = code;
/* Display error message defined
* by the code argument.
*/
if( (code <= 0) || (code >= 7) )
code = 0;
printf( "%s error\n", ermsg[code] );
#endif
/* Return to calling
* program
*/
return( 0 );
}

File diff suppressed because it is too large Load Diff

View File

@@ -0,0 +1,30 @@
acosh.obj,-
asin.obj,-
asinh.obj,-
atan.obj,-
atanh.obj,-
cbrt.obj,-
chbevl.obj,-
const.obj,-
cosh.obj,-
drand.obj,-
exp.obj,-
exp10.obj,-
fabs.obj,-
floor.obj,-
log.obj,-
log10.obj,-
polevl.obj,-
pow.obj,-
powi.obj,-
round.obj,-
sin.obj,-
sinh.obj,-
tan.obj,-
tanh.obj,-
unity.obj,-
sqrt.obj,-
floor.obj,-
polevl.obj,-
mtherr.obj,-
sys$library:vaxcrtl/lib

View File

@@ -0,0 +1,116 @@
; Static Name Aliases
;
TITLE polevl
_TEXT SEGMENT BYTE PUBLIC 'CODE'
_TEXT ENDS
CONST SEGMENT WORD PUBLIC 'CONST'
CONST ENDS
_BSS SEGMENT WORD PUBLIC 'BSS'
_BSS ENDS
_DATA SEGMENT WORD PUBLIC 'DATA'
_DATA ENDS
DGROUP GROUP CONST, _BSS, _DATA
ASSUME CS: _TEXT, DS: DGROUP, SS: DGROUP, ES: DGROUP
PUBLIC _polevl
PUBLIC _p1evl
_DATA SEGMENT
EXTRN __chkstk:NEAR
EXTRN __fac:NEAR
EXTRN __fltused:NEAR
$T20001 DQ 0000000000H ; .0000000000000000
ans DQ 0
ctrlw DW 0
_DATA ENDS
_TEXT SEGMENT
PUBLIC _polevl
_polevl PROC NEAR
push bp
mov bp,sp
mov ax,12
call __chkstk
push si
mov si,[bp+12]
; fstcw ctrlw
; fwait
; mov ax,ctrlw
; or ax,00100h
; mov ctrlw,ax
; fldcw ctrlw
fldz
fwait
mov ax,[bp+14]
inc ax
mov [bp-12],ax
$D15:
fmul QWORD PTR [bp+4]
add si,8
fwait
fadd QWORD PTR [si-8]
fwait
dec WORD PTR [bp-12]
jne $D15
fstp ans
; fstcw ctrlw
; fwait
; mov ax,ctrlw
; and ax,0feffh
; mov ctrlw,ax
; fldcw ctrlw
lea ax, ans
fwait
pop si
mov sp,bp
pop bp
ret
_polevl ENDP
PUBLIC _p1evl
_p1evl PROC NEAR
push bp
mov bp,sp
mov ax,12
call __chkstk
push si
; fstcw ctrlw
; fwait
; mov ax,ctrlw
; or ax,00100h
; mov ctrlw,ax
; fldcw ctrlw
mov si,[bp+12]
fld QWORD PTR [bp+4]
add si,8
fadd QWORD PTR [si-8]
fwait
mov ax,[bp+14]
dec ax
mov [bp-12],ax
$D26:
fmul QWORD PTR [bp+4]
add si,8
fadd QWORD PTR [si-8]
fwait
dec WORD PTR [bp-12]
jne $D26
fstp ans
lea ax, ans
; fstcw ctrlw
; fwait
; mov ax,ctrlw
; and ax,0feffh
; mov ctrlw,ax
; fldcw ctrlw
fwait
pop si
mov sp,bp
pop bp
ret
_p1evl ENDP
_TEXT ENDS
END

View File

@@ -0,0 +1,97 @@
/* polevl.c
* p1evl.c
*
* Evaluate polynomial
*
*
*
* SYNOPSIS:
*
* int N;
* double x, y, coef[N+1], polevl[];
*
* y = polevl( x, coef, N );
*
*
*
* DESCRIPTION:
*
* Evaluates polynomial of degree N:
*
* 2 N
* y = C + C x + C x +...+ C x
* 0 1 2 N
*
* Coefficients are stored in reverse order:
*
* coef[0] = C , ..., coef[N] = C .
* N 0
*
* The function p1evl() assumes that coef[N] = 1.0 and is
* omitted from the array. Its calling arguments are
* otherwise the same as polevl().
*
*
* SPEED:
*
* In the interest of speed, there are no checks for out
* of bounds arithmetic. This routine is used by most of
* the functions in the library. Depending on available
* equipment features, the user may wish to rewrite the
* program in microcode or assembly language.
*
*/
/*
Cephes Math Library Release 2.1: December, 1988
Copyright 1984, 1987, 1988 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/
double polevl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
int i;
double *p;
p = coef;
ans = *p++;
i = N;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}
/* p1evl() */
/* N
* Evaluate polynomial when coefficient of x is 1.0.
* Otherwise same as polevl.
*/
double p1evl( x, coef, N )
double x;
double coef[];
int N;
{
double ans;
double *p;
int i;
p = coef;
ans = x + *p++;
i = N-1;
do
ans = ans * x + *p++;
while( --i );
return( ans );
}

File diff suppressed because it is too large Load Diff

Some files were not shown because too many files have changed in this diff Show More