Steal MSL stuff from sms and fdlibm

This commit is contained in:
Roman Sandu
2025-06-23 20:33:40 +03:00
parent a10e202517
commit bf33740d76
59 changed files with 4782 additions and 220 deletions
@@ -0,0 +1,45 @@
#include "libc/stdlib.h"
void _ExitProcess();
void __destroy_global_chain();
void __kill_critical_regions();
extern void (*_dtors[])(void);
int __aborting = 0;
void (*__atexit_funcs[64])(void);
int __atexit_curr_func = 0;
void (*__stdio_exit)(void) = 0;
void (*__console_exit)(void) = 0;
void exit(int status) {
int i;
void (**dtor)(void);
if (!__aborting) {
__destroy_global_chain();
dtor = _dtors;
while (*dtor != NULL) {
(*dtor)();
dtor++;
}
if (__stdio_exit != NULL) {
__stdio_exit();
__stdio_exit = NULL;
}
}
while (__atexit_curr_func > 0)
__atexit_funcs[--__atexit_curr_func]();
if (__console_exit != NULL) {
__console_exit();
__console_exit = NULL;
}
_ExitProcess();
}
@@ -0,0 +1,131 @@
#include "MSL_C/ansi_files.h"
char stderr_buff[0x100];
char stdout_buff[0x100];
char stdin_buff[0x100];
extern int __close_console(__file_handle file);
extern int __write_console(__file_handle file, unsigned char* buf, size_t* count, __idle_proc idle_fn);
extern int __read_console(__file_handle file, unsigned char* buf, size_t* count, __idle_proc idle_fn);
extern files __files = {
{
0,
0,
1,
1,
2,
0,
0,
0,
0,
0,
0,
0,
0,
{ 0, 0 },
{ 0, 0 },
0,
(unsigned char*)&stdin_buff,
0x100,
(unsigned char*)&stdin_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files._stdout,
},
{
1,
0,
2,
1,
2,
0,
0,
0,
0,
0,
0,
0,
0,
{ 0, 0 },
{ 0, 0 },
0,
(unsigned char*)&stdout_buff,
0x100,
(unsigned char*)&stdout_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files._stderr,
},
{
2,
0,
2,
0,
2,
0,
0,
0,
0,
0,
0,
0,
0,
{ 0, 0 },
{ 0, 0 },
0,
(unsigned char*)&stderr_buff,
0x100,
(unsigned char*)&stderr_buff,
0,
0,
0,
0,
NULL,
__read_console,
__write_console,
__close_console,
NULL,
&__files._sentinel,
},
};
unsigned int __flush_all(void) {
}
void __close_all(void) {
FILE* file = &__files._stdin;
while (file) {
FILE* curr;
if (file->file_mode.file_kind & 0x3)
fclose(file);
curr = file->next;
file = file->next;
if (curr->char_buffer) {
free(curr);
continue;
}
file->file_mode.file_kind = 3;
if (file && file->char_buffer)
curr->next = NULL;
}
}
+343
View File
@@ -0,0 +1,343 @@
#include "MSL_C/ansi_fp.h"
#include "MSL_C/MSL_Common/float.h"
#include "ctype.h"
#include "limits.h"
static int __count_trailing_zerol(unsigned long x) {
int result = 0;
int bits_not_checked = sizeof(unsigned long) * CHAR_BIT;
int n = bits_not_checked / 2;
int mask_size = n;
unsigned long mask = (~0UL) >> (bits_not_checked - n);
while (bits_not_checked) {
if (!(x & mask)) {
result += mask_size;
x >>= mask_size;
bits_not_checked -= mask_size;
} else if (mask == 1) {
break;
}
if (n > 1) {
n /= 2;
}
if (mask > 1) {
mask >>= n;
mask_size -= n;
}
}
return result;
}
static int __count_trailing_zero(double x) {
unsigned long* l = (unsigned long*)&x;
if (l[1] != 0) {
return __count_trailing_zerol(l[1]);
}
return (int)(sizeof(unsigned long) * CHAR_BIT + __count_trailing_zerol(l[0] | 0x00100000));
}
static void __dorounddecup(decimal* d, int digits) {
unsigned char* b = d->sig.text;
unsigned char* i = b + digits - 1;
while (1) {
if (*i < 9) {
*i += 1;
break;
}
if (i == b) {
*i = 1;
d->exp++;
break;
}
*i-- = 0;
}
}
void __ull2dec(decimal* result, unsigned long long val) {
result->sign = 0;
if (val == 0) {
result->exp = 0;
result->sig.length = 1;
result->sig.text[0] = 0;
return;
}
if (val < 0) {
val = -val;
result->sign = 1;
}
result->sig.length = 0;
for (; val != 0; val /= 10) {
result->sig.text[result->sig.length++] = (unsigned char)(val % 10);
}
{
unsigned char* i = result->sig.text;
unsigned char* j = result->sig.text + result->sig.length;
for (; i < --j; ++i) {
unsigned char t = *i;
*i = *j;
*j = t;
}
}
result->exp = result->sig.length - 1;
}
void __timesdec(decimal* result, const decimal* x, const decimal* y) {
unsigned long accumulator = 0;
unsigned char mantissa[SIGDIGLEN * 2];
int i = x->sig.length + y->sig.length - 1;
unsigned char* pDigit;
unsigned char* ip = mantissa + i + 1;
unsigned char* ep = ip;
result->sign = 0;
for (; i > 0; i--) {
int k = y->sig.length - 1;
int j = i - k - 1;
int l;
int t;
const unsigned char* jp;
const unsigned char* kp;
if (j < 0) {
j = 0;
k = i - 1;
}
jp = x->sig.text + j;
kp = y->sig.text + k;
l = k + 1;
t = x->sig.length - j;
if (l > t)
l = t;
for (; l > 0; --l, ++jp, --kp) {
accumulator += *jp * *kp;
}
*--ip = (unsigned char)(accumulator % 10);
accumulator /= 10;
}
result->exp = (short)(x->exp + y->exp);
if (accumulator) {
*--ip = (unsigned char)(accumulator);
result->exp++;
}
for (i = 0; i < SIGDIGLEN && ip < ep; ++i, ++ip) {
result->sig.text[i] = *ip;
}
result->sig.length = (unsigned char)(i);
if (ip < ep && *ip >= 5) {
if (*ip == 5) {
unsigned char* jp = ip + 1;
for (; jp < ep; jp++) {
if (*jp != 0)
goto round;
}
if ((ip[-1] & 1) == 0)
return;
}
round:
__dorounddecup(result, result->sig.length);
}
}
void __str2dec(decimal* d, const char* s, short exp) {
int i;
d->exp = exp;
d->sign = 0;
for (i = 0; i < SIGDIGLEN && *s;) {
d->sig.text[i++] = *s++ - '0';
}
d->sig.length = i;
if (*s != 0) {
if (*s < 5)
return;
if (*s > 5)
goto round;
{
const char* p = s + 1;
for (; *p != 0; p++) {
if (*p != '0')
goto round;
}
if ((d->sig.text[i - 1] & 1) == 0)
return;
}
round:
__dorounddecup(d, d->sig.length);
}
}
void __two_exp(decimal* result, long exp) {
switch (exp) {
case -64:
__str2dec(result, "542101086242752217003726400434970855712890625", -20);
return;
case -53:
__str2dec(result, "11102230246251565404236316680908203125", -16);
return;
case -32:
__str2dec(result, "23283064365386962890625", -10);
return;
case -16:
__str2dec(result, "152587890625", -5);
return;
case -8:
__str2dec(result, "390625", -3);
return;
case -7:
__str2dec(result, "78125", -3);
return;
case -6:
__str2dec(result, "15625", -2);
return;
case -5:
__str2dec(result, "3125", -2);
return;
case -4:
__str2dec(result, "625", -2);
return;
case -3:
__str2dec(result, "125", -1);
return;
case -2:
__str2dec(result, "25", -1);
return;
case -1:
__str2dec(result, "5", -1);
return;
case 0:
__str2dec(result, "1", 0);
return;
case 1:
__str2dec(result, "2", 0);
return;
case 2:
__str2dec(result, "4", 0);
return;
case 3:
__str2dec(result, "8", 0);
return;
case 4:
__str2dec(result, "16", 1);
return;
case 5:
__str2dec(result, "32", 1);
return;
case 6:
__str2dec(result, "64", 1);
return;
case 7:
__str2dec(result, "128", 2);
return;
case 8:
__str2dec(result, "256", 2);
return;
}
{
decimal x2, temp;
__two_exp(&x2, exp / 2);
__timesdec(result, &x2, &x2);
if (exp & 1) {
temp = *result;
if (exp > 0) {
__str2dec(&x2, "2", 0);
} else {
__str2dec(&x2, "5", -1);
}
__timesdec(result, &temp, &x2);
}
}
}
void __num2dec_internal(decimal* d, double x) {
char sign = (char)(signbit(x) != 0);
if (x == 0) {
d->sign = sign;
d->exp = 0;
d->sig.length = 1;
d->sig.text[0] = 0;
return;
}
if (!isfinite(x)) {
d->sign = sign;
d->exp = 0;
d->sig.length = 1;
d->sig.text[0] = fpclassify(x) == 1 ? 'N' : 'I';
return;
}
if (sign != 0) {
x = -x;
}
{
int exp;
double frac = frexp(x, &exp);
long num_bits_extract = DBL_MANT_DIG - __count_trailing_zero(frac);
double integer;
decimal int_d, pow2_d;
__two_exp(&pow2_d, exp - num_bits_extract);
frac = modf(ldexp(frac, num_bits_extract), &integer);
__ull2dec(&int_d, (unsigned long long)integer);
__timesdec(d, &int_d, &pow2_d);
d->sign = sign;
}
}
void __num2dec(const decform* form, double x, decimal* d) {
short digits = form->digits;
int i;
__num2dec_internal(d, x);
if (d->sig.text[0] > 9) {
return;
}
if (digits > SIGDIGLEN) {
digits = SIGDIGLEN;
}
__rounddec(d, digits);
while (d->sig.length < digits) {
d->sig.text[d->sig.length++] = 0;
}
d->exp -= d->sig.length - 1;
for (i = 0; i < d->sig.length; i++) {
d->sig.text[i] += '0';
}
}
@@ -0,0 +1,39 @@
#include "MSL_C/ansi_files.h"
void __prep_buffer(FILE* file) {
file->buffer_ptr = file->buffer;
file->buffer_length = file->buffer_size;
file->buffer_length -= file->position & file->buffer_alignment;
file->buffer_position = file->position;
}
void __convert_from_newlines(unsigned char* p, size_t* n) {
}
int __flush_buffer(FILE* file, size_t* bytes_flushed) {
size_t buffer_len;
int ioresult;
buffer_len = file->buffer_ptr - file->buffer;
if (buffer_len) {
file->buffer_length = buffer_len;
if (!file->file_mode.binary_io)
__convert_from_newlines(file->buffer, &file->buffer_length);
ioresult = (*file->write_fn)(file->handle, file->buffer, &file->buffer_length, file->idle_fn);
if (bytes_flushed)
*bytes_flushed = file->buffer_length;
if (ioresult)
return ioresult;
file->position += file->buffer_length;
}
__prep_buffer(file);
return __no_io_error;
}
+75
View File
@@ -0,0 +1,75 @@
#include "MSL_C/ctype.h"
#define ctrl __control_char
#define motn __motion_char
#define spac __space_char
#define punc __punctuation
#define digi __digit
#define hexd __hex_digit
#define lowc __lower_case
#define uppc __upper_case
#define dhex (hexd | digi)
#define uhex (hexd | uppc)
#define lhex (hexd | lowc)
unsigned char __ctype_map[256] = {
// clang-format off
ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, motn, motn, motn, motn, motn, ctrl, ctrl,
ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl, ctrl,
spac, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc, punc,
dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, dhex, punc, punc, punc, punc, punc, punc,
punc, uhex, uhex, uhex, uhex, uhex, uhex, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc,
uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, uppc, punc, punc, punc, punc, punc,
punc, lhex, lhex, lhex, lhex, lhex, lhex, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc,
lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, lowc, punc, punc, punc, punc, ctrl,
// clang-format on
};
unsigned char __lower_map[256] = {
// clang-format off
0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B, 0x0C, 0x0D, 0x0E, 0x0F,
0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F,
' ', '!', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '.', '/',
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ':', ';', '<', '=', '>', '?',
'@', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', '[', '\\', ']', '^', '_',
'`', 'a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'i', 'j', 'k', 'l', 'm', 'n', 'o',
'p', 'q', 'r', 's', 't', 'u', 'v', 'w', 'x', 'y', 'z', '{', '|', '}', '~', 0x7F,
0x80, 0x81, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x8D, 0x8E, 0x8F,
0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99, 0x9A, 0x9B, 0x9C, 0x9D, 0x9E, 0x9F,
0xA0, 0xA1, 0xA2, 0xA3, 0xA4, 0xA5, 0xA6, 0xA7, 0xA8, 0xA9, 0xAA, 0xAB, 0xAC, 0xAD, 0xAE, 0xAF,
0xB0, 0xB1, 0xB2, 0xB3, 0xB4, 0xB5, 0xB6, 0xB7, 0xB8, 0xB9, 0xBA, 0xBB, 0xBC, 0xBD, 0xBE, 0xBF,
0xC0, 0xC1, 0xC2, 0xC3, 0xC4, 0xC5, 0xC6, 0xC7, 0xC8, 0xC9, 0xCA, 0xCB, 0xCC, 0xCD, 0xCE, 0xCF,
0xD0, 0xD1, 0xD2, 0xD3, 0xD4, 0xD5, 0xD6, 0xD7, 0xD8, 0xD9, 0xDA, 0xDB, 0xDC, 0xDD, 0xDE, 0xDF,
0xE0, 0xE1, 0xE2, 0xE3, 0xE4, 0xE5, 0xE6, 0xE7, 0xE8, 0xE9, 0xEA, 0xEB, 0xEC, 0xED, 0xEE, 0xEF,
0xF0, 0xF1, 0xF2, 0xF3, 0xF4, 0xF5, 0xF6, 0xF7, 0xF8, 0xF9, 0xFA, 0xFB, 0xFC, 0xFD, 0xFE, 0xFF,
// clang-format on
};
unsigned char __upper_map[256] = {
// clang-format off
0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0x08, 0x09, 0x0A, 0x0B, 0x0C, 0x0D, 0x0E, 0x0F,
0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0x18, 0x19, 0x1A, 0x1B, 0x1C, 0x1D, 0x1E, 0x1F,
' ', '!', '"', '#', '$', '%', '&', '\'', '(', ')', '*', '+', ',', '-', '.', '/',
'0', '1', '2', '3', '4', '5', '6', '7', '8', '9', ':', ';', '<', '=', '>', '?',
'@', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '[', '\\', ']', '^', '_',
'`', 'A', 'B', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'J', 'K', 'L', 'M', 'N', 'O',
'P', 'Q', 'R', 'S', 'T', 'U', 'V', 'W', 'X', 'Y', 'Z', '{', '|', '}', '~', 0x7F,
0x80, 0x81, 0x82, 0x83, 0x84, 0x85, 0x86, 0x87, 0x88, 0x89, 0x8A, 0x8B, 0x8C, 0x8D, 0x8E, 0x8F,
0x90, 0x91, 0x92, 0x93, 0x94, 0x95, 0x96, 0x97, 0x98, 0x99, 0x9A, 0x9B, 0x9C, 0x9D, 0x9E, 0x9F,
0xA0, 0xA1, 0xA2, 0xA3, 0xA4, 0xA5, 0xA6, 0xA7, 0xA8, 0xA9, 0xAA, 0xAB, 0xAC, 0xAD, 0xAE, 0xAF,
0xB0, 0xB1, 0xB2, 0xB3, 0xB4, 0xB5, 0xB6, 0xB7, 0xB8, 0xB9, 0xBA, 0xBB, 0xBC, 0xBD, 0xBE, 0xBF,
0xC0, 0xC1, 0xC2, 0xC3, 0xC4, 0xC5, 0xC6, 0xC7, 0xC8, 0xC9, 0xCA, 0xCB, 0xCC, 0xCD, 0xCE, 0xCF,
0xD0, 0xD1, 0xD2, 0xD3, 0xD4, 0xD5, 0xD6, 0xD7, 0xD8, 0xD9, 0xDA, 0xDB, 0xDC, 0xDD, 0xDE, 0xDF,
0xE0, 0xE1, 0xE2, 0xE3, 0xE4, 0xE5, 0xE6, 0xE7, 0xE8, 0xE9, 0xEA, 0xEB, 0xEC, 0xED, 0xEE, 0xEF,
0xF0, 0xF1, 0xF2, 0xF3, 0xF4, 0xF5, 0xF6, 0xF7, 0xF8, 0xF9, 0xFA, 0xFB, 0xFC, 0xFD, 0xFE, 0xFF,
// clang-format on
};
int tolower(int __c) {
if (__c == -1)
return -1;
return __lower_map[__c & 0xff];
}
+114
View File
@@ -0,0 +1,114 @@
#include "MSL_C/direct_io.h"
#include "MSL_C/buffer_io.h"
#include "MSL_C/string.h"
#include "MSL_C/wchar_io.h"
size_t __fwrite(const void* buffer, size_t size, size_t count, FILE* stream);
size_t fwrite(const void* buffer, size_t size, size_t count, FILE* stream) {
__fwrite(buffer, size, count, stream);
}
size_t __fwrite(const void* buffer, size_t size, size_t count, FILE* stream) {
unsigned char* write_ptr;
size_t num_bytes, bytes_to_go, bytes_written;
int ioresult, always_buffer;
if (fwide(stream, 0) == 0)
fwide(stream, -1);
bytes_to_go = size * count;
if (!bytes_to_go || stream->file_state.error || stream->file_mode.file_kind == __closed_file)
return 0;
if (stream->file_mode.file_kind == __console_file)
__stdio_atexit();
always_buffer = !stream->file_mode.binary_io || stream->file_mode.buffer_mode == _IOFBF ||
stream->file_mode.buffer_mode == _IOLBF;
if (stream->file_state.io_state == __neutral && stream->file_mode.io_mode & __write) {
if (stream->file_mode.io_mode & __append) {
if (fseek(stream, 0, 2))
return 0;
}
stream->file_state.io_state = __writing;
__prep_buffer(stream);
}
if (stream->file_state.io_state != __writing) {
set_error(stream);
return 0;
}
write_ptr = (unsigned char*)buffer;
bytes_written = 0;
if (bytes_to_go && (stream->buffer_ptr != stream->buffer || always_buffer)) {
stream->buffer_length = stream->buffer_size - (stream->buffer_ptr - stream->buffer);
do {
unsigned char* newline = NULL;
num_bytes = stream->buffer_length;
if (num_bytes > bytes_to_go)
num_bytes = bytes_to_go;
if (num_bytes) {
memcpy(stream->buffer_ptr, write_ptr, num_bytes);
write_ptr += num_bytes;
bytes_written += num_bytes;
bytes_to_go -= num_bytes;
stream->buffer_ptr += num_bytes;
stream->buffer_length -= num_bytes;
}
if (!stream->buffer_length && (int)stream->file_mode.file_kind == __string_file) {
bytes_written += bytes_to_go;
break;
}
if (!stream->buffer_length || newline != NULL || (stream->file_mode.buffer_mode == _IONBF)) {
ioresult = __flush_buffer(stream, NULL);
if (ioresult) {
set_error(stream);
bytes_to_go = 0;
break;
}
}
} while (bytes_to_go && always_buffer);
}
if (bytes_to_go && !always_buffer) {
unsigned char* save_buffer = stream->buffer;
size_t save_size = stream->buffer_size;
stream->buffer = write_ptr;
stream->buffer_size = bytes_to_go;
stream->buffer_ptr = write_ptr + bytes_to_go;
if (__flush_buffer(stream, &num_bytes) != __no_io_error)
set_error(stream);
bytes_written += num_bytes;
stream->buffer = save_buffer;
stream->buffer_size = save_size;
__prep_buffer(stream);
stream->buffer_length = 0;
}
if (stream->file_mode.buffer_mode != _IOFBF)
stream->buffer_length = 0;
return ((bytes_written + size - 1) / size);
}
+107
View File
@@ -0,0 +1,107 @@
/* @(#)e_acos.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __ieee754_acos(x)
* Method :
* acos(x) = pi/2 - asin(x)
* acos(-x) = pi/2 + asin(x)
* For |x|<=0.5
* acos(x) = pi/2 - (x + x*x^2*R(x^2)) (see asin.c)
* For x>0.5
* acos(x) = pi/2 - (pi/2 - 2asin(sqrt((1-x)/2)))
* = 2asin(sqrt((1-x)/2))
* = 2s + 2s*z*R(z) ...z=(1-x)/2, s=sqrt(z)
* = 2f + (2c + 2s*z*R(z))
* where f=hi part of s, and c = (z-f*f)/(s+f) is the correction term
* for f so that f+c ~ sqrt(z).
* For x<-0.5
* acos(x) = pi - 2asin(sqrt((1-|x|)/2))
* = pi - 0.5*(s+s*z*R(z)), where z=(1-|x|)/2,s=sqrt(z)
*
* Special cases:
* if x is NaN, return x itself;
* if |x|>1, return NaN with invalid signal.
*
* Function needed: sqrt
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double
#else
static double
#endif
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
pi = 3.14159265358979311600e+00, /* 0x400921FB, 0x54442D18 */
pio2_hi = 1.57079632679489655800e+00, /* 0x3FF921FB, 0x54442D18 */
pio2_lo = 6.12323399573676603587e-17, /* 0x3C91A626, 0x33145C07 */
pS0 = 1.66666666666666657415e-01, /* 0x3FC55555, 0x55555555 */
pS1 = -3.25565818622400915405e-01, /* 0xBFD4D612, 0x03EB6F7D */
pS2 = 2.01212532134862925881e-01, /* 0x3FC9C155, 0x0E884455 */
pS3 = -4.00555345006794114027e-02, /* 0xBFA48228, 0xB5688F3B */
pS4 = 7.91534994289814532176e-04, /* 0x3F49EFE0, 0x7501B288 */
pS5 = 3.47933107596021167570e-05, /* 0x3F023DE1, 0x0DFDF709 */
qS1 = -2.40339491173441421878e+00, /* 0xC0033A27, 0x1C8A2D4B */
qS2 = 2.02094576023350569471e+00, /* 0x40002AE5, 0x9C598AC8 */
qS3 = -6.88283971605453293030e-01, /* 0xBFE6066C, 0x1B8D0159 */
qS4 = 7.70381505559019352791e-02; /* 0x3FB3B8C5, 0xB12E9282 */
#ifdef __STDC__
double __ieee754_acos(double x)
#else
double __ieee754_acos(x)
double x;
#endif
{
double z, p, q, r, w, s, c, df;
int hx, ix;
hx = __HI(x);
ix = hx & 0x7fffffff;
if (ix >= 0x3ff00000) { /* |x| >= 1 */
if (((ix - 0x3ff00000) | __LO(x)) == 0) { /* |x|==1 */
if (hx > 0)
return 0.0; /* acos(1) = 0 */
else
return pi + 2.0 * pio2_lo; /* acos(-1)= pi */
}
return NAN; /* acos(|x|>1) is NaN */
}
if (ix < 0x3fe00000) { /* |x| < 0.5 */
if (ix <= 0x3c600000)
return pio2_hi + pio2_lo; /*if|x|<2**-57*/
z = x * x;
p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
r = p / q;
return pio2_hi - (x - (pio2_lo - x * r));
} else if (hx < 0) { /* x < -0.5 */
z = (one + x) * 0.5;
p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
s = sqrt(z);
r = p / q;
w = r * s - pio2_lo;
return pi - 2.0 * (s + w);
} else { /* x > 0.5 */
z = (one - x) * 0.5;
s = sqrt(z);
df = s;
__LO(df) = 0;
c = (z - df * df) / (s + df);
p = z * (pS0 + z * (pS1 + z * (pS2 + z * (pS3 + z * (pS4 + z * pS5)))));
q = one + z * (qS1 + z * (qS2 + z * (qS3 + z * qS4)));
r = p / q;
w = r * s + c;
return 2.0 * (df + w);
}
}
+143
View File
@@ -0,0 +1,143 @@
/* @(#)e_atan2.c 1.3 95/01/18 */
/**
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
/* __ieee754_atan2(y,x)
* Method :
* 1. Reduce y to positive by atan2(y,x)=-atan2(-y,x).
* 2. Reduce x to positive by (if x and y are unexceptional):
* ARG (x+iy) = arctan(y/x) ... if x > 0,
* ARG (x+iy) = pi - arctan[y/(-x)] ... if x < 0,
*
* Special cases:
*
* ATAN2((anything), NaN ) is NaN;
* ATAN2(NAN , (anything) ) is NaN;
* ATAN2(+-0, +(anything but NaN)) is +-0 ;
* ATAN2(+-0, -(anything but NaN)) is +-pi ;
* ATAN2(+-(anything but 0 and NaN), 0) is +-pi/2;
* ATAN2(+-(anything but INF and NaN), +INF) is +-0 ;
* ATAN2(+-(anything but INF and NaN), -INF) is +-pi;
* ATAN2(+-INF,+INF ) is +-pi/4 ;
* ATAN2(+-INF,-INF ) is +-3pi/4;
* ATAN2(+-INF, (anything but,0,NaN, and INF)) is +-pi/2;
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double
#else
static double
#endif
tiny = 1.0e-300,
zero = 0.0, pi_o_4 = 7.8539816339744827900E-01, /* 0x3FE921FB, 0x54442D18 */
pi_o_2 = 1.5707963267948965580E+00, /* 0x3FF921FB, 0x54442D18 */
pi = 3.1415926535897931160E+00, /* 0x400921FB, 0x54442D18 */
pi_lo = 1.2246467991473531772E-16; /* 0x3CA1A626, 0x33145C07 */
#ifdef __STDC__
double __ieee754_atan2(double y, double x)
#else
double __ieee754_atan2(y, x)
double y, x;
#endif
{
double z;
int k, m, hx, hy, ix, iy;
unsigned lx, ly;
hx = __HI(x);
ix = hx & 0x7fffffff;
lx = __LO(x);
hy = __HI(y);
iy = hy & 0x7fffffff;
ly = __LO(y);
if (((ix | ((lx | -lx) >> 31)) > 0x7ff00000) || ((iy | ((ly | -ly) >> 31)) > 0x7ff00000)) /* x or y is NaN */
return x + y;
if ((hx - 0x3ff00000 | lx) == 0)
return atan(y); /* x=1.0 */
m = ((hy >> 31) & 1) | ((hx >> 30) & 2); /* 2*sign(x)+sign(y) */
/* when y = 0 */
if ((iy | ly) == 0) {
switch (m) {
case 0:
case 1:
return y; /* atan(+-0,+anything)=+-0 */
case 2:
return pi + tiny; /* atan(+0,-anything) = pi */
case 3:
return -pi - tiny; /* atan(-0,-anything) =-pi */
}
}
/* when x = 0 */
if ((ix | lx) == 0)
return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
/* when x is INF */
if (ix == 0x7ff00000) {
if (iy == 0x7ff00000) {
switch (m) {
case 0:
return pi_o_4 + tiny; /* atan(+INF,+INF) */
case 1:
return -pi_o_4 - tiny; /* atan(-INF,+INF) */
case 2:
return 3.0 * pi_o_4 + tiny; /*atan(+INF,-INF)*/
case 3:
return -3.0 * pi_o_4 - tiny; /*atan(-INF,-INF)*/
}
} else {
switch (m) {
case 0:
return zero; /* atan(+...,+INF) */
case 1:
return -zero; /* atan(-...,+INF) */
case 2:
return pi + tiny; /* atan(+...,-INF) */
case 3:
return -pi - tiny; /* atan(-...,-INF) */
}
}
}
/* when y is INF */
if (iy == 0x7ff00000)
return (hy < 0) ? -pi_o_2 - tiny : pi_o_2 + tiny;
/* compute y/x */
k = (iy - ix) >> 20;
if (k > 60)
z = pi_o_2 + 0.5 * pi_lo; /* |y/x| > 2**60 */
else if (hx < 0 && k < -60)
z = 0.0; /* |y|/x < -2**60 */
else
z = atan(fabs(y / x)); /* safe to do y/x */
switch (m) {
case 0:
return z; /* atan(+,+) */
case 1:
__HI(z) ^= 0x80000000;
return z; /* atan(-,+) */
case 2:
return pi - (z - pi_lo); /* atan(+,-) */
default: /* case 3 */
return (z - pi_lo) - pi; /* atan(-,-) */
}
}
@@ -0,0 +1,181 @@
/* @(#)e_rem_pio2.c 1.4 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
/* __ieee754_rem_pio2(x,y)
*
* return the remainder of x rem pi/2 in y[0]+y[1]
* use __kernel_rem_pio2()
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
/*
* Table of constants for 2/pi, 396 Hex digits (476 decimal) of 2/pi
*/
#ifdef __STDC__
static const int two_over_pi[] = {
#else
static int two_over_pi[] = {
#endif
0xA2F983, 0x6E4E44, 0x1529FC, 0x2757D1, 0xF534DD, 0xC0DB62, 0x95993C, 0x439041, 0xFE5163, 0xABDEBB, 0xC561B7,
0x246E3A, 0x424DD2, 0xE00649, 0x2EEA09, 0xD1921C, 0xFE1DEB, 0x1CB129, 0xA73EE8, 0x8235F5, 0x2EBB44, 0x84E99C,
0x7026B4, 0x5F7E41, 0x3991D6, 0x398353, 0x39F49C, 0x845F8B, 0xBDF928, 0x3B1FF8, 0x97FFDE, 0x05980F, 0xEF2F11,
0x8B5A0A, 0x6D1F6D, 0x367ECF, 0x27CB09, 0xB74F46, 0x3F669E, 0x5FEA2D, 0x7527BA, 0xC7EBE5, 0xF17B3D, 0x0739F7,
0x8A5292, 0xEA6BFB, 0x5FB11F, 0x8D5D08, 0x560330, 0x46FC7B, 0x6BABF0, 0xCFBC20, 0x9AF436, 0x1DA9E3, 0x91615E,
0xE61B08, 0x659985, 0x5F14A0, 0x68408D, 0xFFD880, 0x4D7327, 0x310606, 0x1556CA, 0x73A8C9, 0x60E27B, 0xC08C6B,
};
#ifdef __STDC__
static const int npio2_hw[] = {
#else
static int npio2_hw[] = {
#endif
0x3FF921FB, 0x400921FB, 0x4012D97C, 0x401921FB, 0x401F6A7A, 0x4022D97C, 0x4025FDBB, 0x402921FB,
0x402C463A, 0x402F6A7A, 0x4031475C, 0x4032D97C, 0x40346B9C, 0x4035FDBB, 0x40378FDB, 0x403921FB,
0x403AB41B, 0x403C463A, 0x403DD85A, 0x403F6A7A, 0x40407E4C, 0x4041475C, 0x4042106C, 0x4042D97C,
0x4043A28C, 0x40446B9C, 0x404534AC, 0x4045FDBB, 0x4046C6CB, 0x40478FDB, 0x404858EB, 0x404921FB,
};
/*
* invpio2: 53 bits of 2/pi
* pio2_1: first 33 bit of pi/2
* pio2_1t: pi/2 - pio2_1
* pio2_2: second 33 bit of pi/2
* pio2_2t: pi/2 - (pio2_1+pio2_2)
* pio2_3: third 33 bit of pi/2
* pio2_3t: pi/2 - (pio2_1+pio2_2+pio2_3)
*/
#ifdef __STDC__
static const double
#else
static double
#endif
zero = 0.00000000000000000000e+00, /* 0x00000000, 0x00000000 */
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
invpio2 = 6.36619772367581382433e-01, /* 0x3FE45F30, 0x6DC9C883 */
pio2_1 = 1.57079632673412561417e+00, /* 0x3FF921FB, 0x54400000 */
pio2_1t = 6.07710050650619224932e-11, /* 0x3DD0B461, 0x1A626331 */
pio2_2 = 6.07710050630396597660e-11, /* 0x3DD0B461, 0x1A600000 */
pio2_2t = 2.02226624879595063154e-21, /* 0x3BA3198A, 0x2E037073 */
pio2_3 = 2.02226624871116645580e-21, /* 0x3BA3198A, 0x2E000000 */
pio2_3t = 8.47842766036889956997e-32; /* 0x397B839A, 0x252049C1 */
#ifdef __STDC__
int __ieee754_rem_pio2(double x, double* y)
#else
int __ieee754_rem_pio2(x, y)
double x, y[];
#endif
{
double z, w, t, r, fn;
double tx[3];
int e0, i, j, nx, n, ix, hx;
hx = __HI(x); /* high word of x */
ix = hx & 0x7fffffff;
if (ix <= 0x3fe921fb) /* |x| ~<= pi/4 , no need for reduction */
{
y[0] = x;
y[1] = 0;
return 0;
}
if (ix < 0x4002d97c) { /* |x| < 3pi/4, special case with n=+-1 */
if (hx > 0) {
z = x - pio2_1;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z - pio2_1t;
y[1] = (z - y[0]) - pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z -= pio2_2;
y[0] = z - pio2_2t;
y[1] = (z - y[0]) - pio2_2t;
}
return 1;
} else { /* negative x */
z = x + pio2_1;
if (ix != 0x3ff921fb) { /* 33+53 bit pi is good enough */
y[0] = z + pio2_1t;
y[1] = (z - y[0]) + pio2_1t;
} else { /* near pi/2, use 33+33+53 bit pi */
z += pio2_2;
y[0] = z + pio2_2t;
y[1] = (z - y[0]) + pio2_2t;
}
return -1;
}
}
if (ix <= 0x413921fb) { /* |x| ~<= 2^19*(pi/2), medium size */
t = fabs(x);
n = (int)(t * invpio2 + half);
fn = (double)n;
r = t - fn * pio2_1;
w = fn * pio2_1t; /* 1st round good to 85 bit */
if (n < 32 && ix != npio2_hw[n - 1]) {
y[0] = r - w; /* quick check no cancellation */
} else {
j = ix >> 20;
y[0] = r - w;
i = j - (((__HI(y[0])) >> 20) & 0x7ff);
if (i > 16) { /* 2nd iteration needed, good to 118 */
t = r;
w = fn * pio2_2;
r = t - w;
w = fn * pio2_2t - ((t - r) - w);
y[0] = r - w;
i = j - (((__HI(y[0])) >> 20) & 0x7ff);
if (i > 49) { /* 3rd iteration need, 151 bits acc */
t = r; /* will cover all possible cases */
w = fn * pio2_3;
r = t - w;
w = fn * pio2_3t - ((t - r) - w);
y[0] = r - w;
}
}
}
y[1] = (r - y[0]) - w;
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
} else
return n;
}
/*
* all other (large) arguments
*/
if (ix >= 0x7ff00000) { /* x is inf or NaN */
y[0] = y[1] = x - x;
return 0;
}
/* set z = scalbn(|x|,ilogb(x)-23) */
__LO(z) = __LO(x);
e0 = (ix >> 20) - 1046; /* e0 = ilogb(z)-23; */
__HI(z) = ix - (e0 << 20);
for (i = 0; i < 2; i++) {
tx[i] = (double)((int)(z));
z = (z - tx[i]) * two24;
}
tx[2] = z;
nx = 3;
while (tx[nx - 1] == zero)
nx--; /* skip zero term */
n = __kernel_rem_pio2(tx, y, e0, nx, 2, two_over_pi);
if (hx < 0) {
y[0] = -y[0];
y[1] = -y[1];
return -n;
}
return n;
}
+3
View File
@@ -0,0 +1,3 @@
#include "errno.h"
int errno;
+2
View File
@@ -0,0 +1,2 @@
unsigned long __float_nan[] = { 0x7FFFFFFF };
unsigned long __float_huge[] = { 0x7F800000 };
+92
View File
@@ -0,0 +1,92 @@
/* @(#)k_cos.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __kernel_cos( x, y )
* kernel cos function on [-pi/4, pi/4], pi/4 ~ 0.785398164
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
*
* Algorithm
* 1. Since cos(-x) = cos(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return 1 with inexact if x!=0.
* 3. cos(x) is approximated by a polynomial of degree 14 on
* [0,pi/4]
* 4 14
* cos(x) ~ 1 - x*x/2 + C1*x + ... + C6*x
* where the remez error is
*
* | 2 4 6 8 10 12 14 | -58
* |cos(x)-(1-.5*x +C1*x +C2*x +C3*x +C4*x +C5*x +C6*x )| <= 2
* | |
*
* 4 6 8 10 12 14
* 4. let r = C1*x +C2*x +C3*x +C4*x +C5*x +C6*x , then
* cos(x) = 1 - x*x/2 + r
* since cos(x+y) ~ cos(x) - sin(x)*y
* ~ cos(x) - x*y,
* a correction term is necessary in cos(x) and hence
* cos(x+y) = 1 - (x*x/2 - (r - x*y))
* For better accuracy when x > 0.3, let qx = |x|/4 with
* the last 32 bits mask off, and if x > 0.78125, let qx = 0.28125.
* Then
* cos(x+y) = (1-qx) - ((x*x/2-qx) - (r-x*y)).
* Note that 1-qx and (x*x/2-qx) is EXACT here, and the
* magnitude of the latter is at least a quarter of x*x/2,
* thus, reducing the rounding error in the subtraction.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double
#else
static double
#endif
one = 1.00000000000000000000e+00, /* 0x3FF00000, 0x00000000 */
C1 = 4.16666666666666019037e-02, /* 0x3FA55555, 0x5555554C */
C2 = -1.38888888888741095749e-03, /* 0xBF56C16C, 0x16C15177 */
C3 = 2.48015872894767294178e-05, /* 0x3EFA01A0, 0x19CB1590 */
C4 = -2.75573143513906633035e-07, /* 0xBE927E4F, 0x809C52AD */
C5 = 2.08757232129817482790e-09, /* 0x3E21EE9E, 0xBDB4B1C4 */
C6 = -1.13596475577881948265e-11; /* 0xBDA8FAE9, 0xBE8838D4 */
#ifdef __STDC__
double __kernel_cos(double x, double y)
#else
double __kernel_cos(x, y)
double x, y;
#endif
{
double a, hz, z, r, qx;
int ix;
ix = __HI(x) & 0x7fffffff; /* ix = |x|'s high word*/
if (ix < 0x3e400000) { /* if x < 2**27 */
if (((int)x) == 0)
return one; /* generate inexact */
}
z = x * x;
r = z * (C1 + z * (C2 + z * (C3 + z * (C4 + z * (C5 + z * C6)))));
if (ix < 0x3FD33333) /* if |x| < 0.3 */
return one - (0.5 * z - (z * r - x * y));
else {
if (ix > 0x3fe90000) { /* x > 0.78125 */
qx = 0.28125;
} else {
__HI(qx) = ix - 0x00200000; /* x/4 */
__LO(qx) = 0;
}
hz = 0.5 * z - qx;
a = one - qx;
return a - (hz - (z * r - x * y));
}
}
@@ -0,0 +1,348 @@
/* @(#)k_rem_pio2.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* __kernel_rem_pio2(x,y,e0,nx,prec,ipio2)
* double x[],y[]; int e0,nx,prec; int ipio2[];
*
* __kernel_rem_pio2 return the last three digits of N with
* y = x - N*pi/2
* so that |y| < pi/2.
*
* The method is to compute the integer (mod 8) and fraction parts of
* (2/pi)*x without doing the full multiplication. In general we
* skip the part of the product that are known to be a huge integer (
* more accurately, = 0 mod 8 ). Thus the number of operations are
* independent of the exponent of the input.
*
* (2/pi) is represented by an array of 24-bit integers in ipio2[].
*
* Input parameters:
* x[] The input value (must be positive) is broken into nx
* pieces of 24-bit integers in double precision format.
* x[i] will be the i-th 24 bit of x. The scaled exponent
* of x[0] is given in input parameter e0 (i.e., x[0]*2^e0
* match x's up to 24 bits.
*
* Example of breaking a double positive z into x[0]+x[1]+x[2]:
* e0 = ilogb(z)-23
* z = scalbn(z,-e0)
* for i = 0,1,2
* x[i] = floor(z)
* z = (z-x[i])*2**24
*
*
* y[] ouput result in an array of double precision numbers.
* The dimension of y[] is:
* 24-bit precision 1
* 53-bit precision 2
* 64-bit precision 2
* 113-bit precision 3
* The actual value is the sum of them. Thus for 113-bit
* precison, one may have to do something like:
*
* long double t,w,r_head, r_tail;
* t = (long double)y[2] + (long double)y[1];
* w = (long double)y[0];
* r_head = t+w;
* r_tail = w - (r_head - t);
*
* e0 The exponent of x[0]
*
* nx dimension of x[]
*
* prec an integer indicating the precision:
* 0 24 bits (single)
* 1 53 bits (double)
* 2 64 bits (extended)
* 3 113 bits (quad)
*
* ipio2[]
* integer array, contains the (24*i)-th to (24*i+23)-th
* bit of 2/pi after binary point. The corresponding
* floating value is
*
* ipio2[i] * 2^(-24(i+1)).
*
* External function:
* double scalbn(), floor();
*
*
* Here is the description of some local variables:
*
* jk jk+1 is the initial number of terms of ipio2[] needed
* in the computation. The recommended value is 2,3,4,
* 6 for single, double, extended,and quad.
*
* jz local integer variable indicating the number of
* terms of ipio2[] used.
*
* jx nx - 1
*
* jv index for pointing to the suitable ipio2[] for the
* computation. In general, we want
* ( 2^e0*x[0] * ipio2[jv-1]*2^(-24jv) )/8
* is an integer. Thus
* e0-3-24*jv >= 0 or (e0-3)/24 >= jv
* Hence jv = max(0,(e0-3)/24).
*
* jp jp+1 is the number of terms in PIo2[] needed, jp = jk.
*
* q[] double array with integral value, representing the
* 24-bits chunk of the product of x and 2/pi.
*
* q0 the corresponding exponent of q[0]. Note that the
* exponent for q[i] would be q0-24*i.
*
* PIo2[] double precision array, obtained by cutting pi/2
* into 24 bits chunks.
*
* f[] ipio2[] in floating point
*
* iq[] integer array by breaking up q[] in 24-bits chunk.
*
* fq[] final product of x*(2/pi) in fq[0],..,fq[jk]
*
* ih integer. If >0 it indicates q[] is >= 0.5, hence
* it also indicates the *sign* of the result.
*
*/
/*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const int init_jk[] = { 2, 3, 4, 6 }; /* initial value for jk */
#else
static int init_jk[] = { 2, 3, 4, 6 };
#endif
#ifdef __STDC__
static const double PIo2[] = {
#else
static double PIo2[] = {
#endif
1.57079625129699707031e+00, /* 0x3FF921FB, 0x40000000 */
7.54978941586159635335e-08, /* 0x3E74442D, 0x00000000 */
5.39030252995776476554e-15, /* 0x3CF84698, 0x80000000 */
3.28200341580791294123e-22, /* 0x3B78CC51, 0x60000000 */
1.27065575308067607349e-29, /* 0x39F01B83, 0x80000000 */
1.22933308981111328932e-36, /* 0x387A2520, 0x40000000 */
2.73370053816464559624e-44, /* 0x36E38222, 0x80000000 */
2.16741683877804819444e-51, /* 0x3569F31D, 0x00000000 */
};
#ifdef __STDC__
static const double
#else
static double
#endif
zero = 0.0,
one = 1.0, two24 = 1.67772160000000000000e+07, /* 0x41700000, 0x00000000 */
twon24 = 5.96046447753906250000e-08; /* 0x3E700000, 0x00000000 */
#ifdef __STDC__
int __kernel_rem_pio2(double* x, double* y, int e0, int nx, int prec, const int* ipio2)
#else
int __kernel_rem_pio2(x, y, e0, nx, prec, ipio2)
double x[], y[];
int e0, nx, prec;
int ipio2[];
#endif
{
int jz, jx, jv, jp, jk, carry, n, iq[20], i, j, k, m, q0, ih;
double z, fw, f[20], fq[20], q[20];
/* initialize jk*/
jk = init_jk[prec];
jp = jk;
/* determine jx,jv,q0, note that 3>q0 */
jx = nx - 1;
jv = (e0 - 3) / 24;
if (jv < 0)
jv = 0;
q0 = e0 - 24 * (jv + 1);
/* set up f[0] to f[jx+jk] where f[jx+jk] = ipio2[jv+jk] */
j = jv - jx;
m = jx + jk;
for (i = 0; i <= m; i++, j++)
f[i] = (j < 0) ? zero : (double)ipio2[j];
/* compute q[0],q[1],...q[jk] */
for (i = 0; i <= jk; i++) {
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
jz = jk;
recompute:
/* distill q[] into iq[] reversingly */
for (i = 0, j = jz, z = q[jz]; j > 0; i++, j--) {
fw = (double)((int)(twon24 * z));
iq[i] = (int)(z - two24 * fw);
z = q[j - 1] + fw;
}
/* compute n */
z = ldexp(z, q0); /* actual value of z */
z -= 8.0 * floor(z * 0.125); /* trim off integer >= 8 */
n = (int)z;
z -= (double)n;
ih = 0;
if (q0 > 0) { /* need iq[jz-1] to determine n */
i = (iq[jz - 1] >> (24 - q0));
n += i;
iq[jz - 1] -= i << (24 - q0);
ih = iq[jz - 1] >> (23 - q0);
} else if (q0 == 0)
ih = iq[jz - 1] >> 23;
else if (z >= 0.5)
ih = 2;
if (ih > 0) { /* q > 0.5 */
n += 1;
carry = 0;
for (i = 0; i < jz; i++) { /* compute 1-q */
j = iq[i];
if (carry == 0) {
if (j != 0) {
carry = 1;
iq[i] = 0x1000000 - j;
}
} else
iq[i] = 0xffffff - j;
}
if (q0 > 0) { /* rare case: chance is 1 in 12 */
switch (q0) {
case 1:
iq[jz - 1] &= 0x7fffff;
break;
case 2:
iq[jz - 1] &= 0x3fffff;
break;
}
}
if (ih == 2) {
z = one - z;
if (carry != 0)
z -= ldexp(one, q0);
}
}
/* check if recomputation is needed */
if (z == zero) {
j = 0;
for (i = jz - 1; i >= jk; i--)
j |= iq[i];
if (j == 0) { /* need recomputation */
for (k = 1; iq[jk - k] == 0; k++)
; /* k = no. of terms needed */
for (i = jz + 1; i <= jz + k; i++) { /* add q[jz+1] to q[jz+k] */
f[jx + i] = (double)ipio2[jv + i];
for (j = 0, fw = 0.0; j <= jx; j++)
fw += x[j] * f[jx + i - j];
q[i] = fw;
}
jz += k;
goto recompute;
}
}
/* chop off zero terms */
if (z == 0.0) {
jz -= 1;
q0 -= 24;
while (iq[jz] == 0) {
jz--;
q0 -= 24;
}
} else { /* break z into 24-bit if necessary */
z = ldexp(z, -q0);
if (z >= two24) {
fw = (double)((int)(twon24 * z));
iq[jz] = (int)(z - two24 * fw);
jz += 1;
q0 += 24;
iq[jz] = (int)fw;
} else
iq[jz] = (int)z;
}
/* convert integer "bit" chunk to floating-point value */
fw = ldexp(one, q0);
for (i = jz; i >= 0; i--) {
q[i] = fw * (double)iq[i];
fw *= twon24;
}
/* compute PIo2[0,...,jp]*q[jz,...,0] */
for (i = jz; i >= 0; i--) {
for (fw = 0.0, k = 0; k <= jp && k <= jz - i; k++)
fw += PIo2[k] * q[i + k];
fq[jz - i] = fw;
}
/* compress fq[] into y[] */
switch (prec) {
case 0:
fw = 0.0;
for (i = jz; i >= 0; i--)
fw += fq[i];
y[0] = (ih == 0) ? fw : -fw;
break;
case 1:
case 2:
fw = 0.0;
for (i = jz; i >= 0; i--)
fw += fq[i];
y[0] = (ih == 0) ? fw : -fw;
fw = fq[0] - fw;
for (i = 1; i <= jz; i++)
fw += fq[i];
y[1] = (ih == 0) ? fw : -fw;
break;
case 3: /* painful */
for (i = jz; i > 0; i--) {
fw = fq[i - 1] + fq[i];
fq[i] += fq[i - 1] - fw;
fq[i - 1] = fw;
}
for (i = jz; i > 1; i--) {
fw = fq[i - 1] + fq[i];
fq[i] += fq[i - 1] - fw;
fq[i - 1] = fw;
}
for (fw = 0.0, i = jz; i >= 2; i--)
fw += fq[i];
if (ih == 0) {
y[0] = fq[0];
y[1] = fq[1];
y[2] = fw;
} else {
y[0] = -fq[0];
y[1] = -fq[1];
y[2] = -fw;
}
}
return n & 7;
}
+79
View File
@@ -0,0 +1,79 @@
/* @(#)k_sin.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* __kernel_sin( x, y, iy)
* kernel sin function on [-pi/4, pi/4], pi/4 ~ 0.7854
* Input x is assumed to be bounded by ~pi/4 in magnitude.
* Input y is the tail of x.
* Input iy indicates whether y is 0. (if iy=0, y assume to be 0).
*
* Algorithm
* 1. Since sin(-x) = -sin(x), we need only to consider positive x.
* 2. if x < 2^-27 (hx<0x3e400000 0), return x with inexact if x!=0.
* 3. sin(x) is approximated by a polynomial of degree 13 on
* [0,pi/4]
* 3 13
* sin(x) ~ x + S1*x + ... + S6*x
* where
*
* |sin(x) 2 4 6 8 10 12 | -58
* |----- - (1+S1*x +S2*x +S3*x +S4*x +S5*x +S6*x )| <= 2
* | x |
*
* 4. sin(x+y) = sin(x) + sin'(x')*y
* ~ sin(x) + (1-x*x/2)*y
* For better accuracy, let
* 3 2 2 2 2
* r = x *(S2+x *(S3+x *(S4+x *(S5+x *S6))))
* then 3 2
* sin(x) = x + (S1*x + (x *(r-y/2)+y))
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double
#else
static double
#endif
half = 5.00000000000000000000e-01, /* 0x3FE00000, 0x00000000 */
S1 = -1.66666666666666324348e-01, /* 0xBFC55555, 0x55555549 */
S2 = 8.33333333332248946124e-03, /* 0x3F811111, 0x1110F8A6 */
S3 = -1.98412698298579493134e-04, /* 0xBF2A01A0, 0x19C161D5 */
S4 = 2.75573137070700676789e-06, /* 0x3EC71DE3, 0x57B1FE7D */
S5 = -2.50507602534068634195e-08, /* 0xBE5AE5E6, 0x8A2B9CEB */
S6 = 1.58969099521155010221e-10; /* 0x3DE5D93A, 0x5ACFD57C */
#ifdef __STDC__
double __kernel_sin(double x, double y, int iy)
#else
double __kernel_sin(x, y, iy)
double x, y;
int iy; /* iy=0 if y is zero */
#endif
{
double z, r, v;
int ix;
ix = __HI(x) & 0x7fffffff; /* high word of x */
if (ix < 0x3e400000) /* |x| < 2**-27 */
{
if ((int)x == 0)
return x;
} /* generate inexact */
z = x * x;
v = z * x;
r = S2 + z * (S3 + z * (S4 + z * (S5 + z * S6)));
if (iy == 0)
return x + v * (S1 + z * r);
else
return x - ((z * (half * y - v * r) - y) - v * S1);
}
@@ -0,0 +1,63 @@
#include "MSL_C/mbstring.h"
#include "MSL_C/string.h"
#pragma dont_inline on
size_t unicode_to_UTF8(char* r3, wchar_t r4) {
const char lut[] = { 0, 0, 0xC0, 0xE0 };
int r5;
char* r6;
if (!r3)
return 0;
if (r4 < 0x80)
r5 = 1;
else if (r4 < 0x800)
r5 = 2;
else
r5 = 3;
r6 = r3 + r5;
switch (r5) {
case 3:
*--r6 = (r4 & 0x3F) | 0x80;
r4 >>= 6;
case 2:
*--r6 = (r4 & 0x3F) | 0x80;
r4 >>= 6;
case 1:
*--r6 = r4 | lut[r5];
}
return r5;
}
#pragma dont_inline off
size_t wcstombs(char* s, const wchar_t* pwcs, size_t n) {
char decoded[4];
wchar_t w;
size_t num;
int i = 0;
const wchar_t* p;
if (!s || !pwcs)
return 0;
p = pwcs;
while (i <= n) {
w = *p;
if (!w) {
s[i] = 0;
break;
}
++p;
num = unicode_to_UTF8(decoded, w);
if (i + num > n)
break;
strncpy(s + i, decoded, num);
i += num;
}
return i;
}
+86
View File
@@ -0,0 +1,86 @@
#include "string.h"
void __copy_longs_rev_unaligned(void* dst, const void* src, size_t n);
void __copy_longs_unaligned(void* dst, const void* src, size_t n);
void __copy_longs_rev_aligned(void* dst, const void* src, size_t n);
void __copy_longs_aligned(void* dst, const void* src, size_t n);
void* memmove(void* dst, const void* src, size_t n) {
unsigned char* csrc;
unsigned char* cdst;
int reverse = (unsigned int)src < (unsigned int)dst;
if (n >= 32) {
if (((unsigned int)dst ^ (unsigned int)src) & 3) {
if (!reverse) {
__copy_longs_unaligned(dst, src, n);
} else {
__copy_longs_rev_unaligned(dst, src, n);
}
} else {
if (!reverse) {
__copy_longs_aligned(dst, src, n);
} else {
__copy_longs_rev_aligned(dst, src, n);
}
}
return dst;
} else {
if (!reverse) {
csrc = ((unsigned char*)src) - 1;
cdst = ((unsigned char*)dst) - 1;
n++;
while (--n > 0) {
*++cdst = *++csrc;
}
} else {
csrc = (unsigned char*)src + n;
cdst = (unsigned char*)dst + n;
n++;
while (--n > 0) {
*--cdst = *--csrc;
}
}
}
return dst;
}
void* memchr(const void* ptr, int ch, size_t count) {
const unsigned char* p;
unsigned long v = (ch & 0xff);
for (p = (unsigned char*)ptr - 1, count++; --count;)
if ((*++p & 0xff) == v)
return (void*)p;
return NULL;
}
void* __memrchr(const void* ptr, int ch, size_t count) {
const unsigned char* p;
unsigned long v = (ch & 0xff);
for (p = (unsigned char*)ptr + count, count++; --count;)
if ((*--p & 0xff) == v)
return (void*)p;
return NULL;
}
int memcmp(const void* lhs, const void* rhs, size_t count) {
const unsigned char* p1;
const unsigned char* p2;
for (p1 = (const unsigned char*)lhs - 1, p2 = (const unsigned char*)rhs - 1, count++; --count;)
if (*++p1 != *++p2)
return ((*p1 < *p2) ? -1 : +1);
return 0;
}
+217
View File
@@ -0,0 +1,217 @@
#include "libc/stddef.h"
#define cps ((unsigned char*)src)
#define cpd ((unsigned char*)dst)
#define lps ((unsigned long*)src)
#define lpd ((unsigned long*)dst)
#define deref_auto_inc(p) *++(p)
void __copy_longs_aligned(void* dst, const void* src, size_t n) {
unsigned long i;
i = (-(unsigned long)dst) & 3;
cps = ((unsigned char*)src) - 1;
cpd = ((unsigned char*)dst) - 1;
if (i) {
n -= i;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--i);
}
lps = ((unsigned long*)(cps + 1)) - 1;
lpd = ((unsigned long*)(cpd + 1)) - 1;
i = n >> 5;
if (i)
do {
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
deref_auto_inc(lpd) = deref_auto_inc(lps);
} while (--i);
i = (n & 31) >> 2;
if (i)
do
deref_auto_inc(lpd) = deref_auto_inc(lps);
while (--i);
cps = ((unsigned char*)(lps + 1)) - 1;
cpd = ((unsigned char*)(lpd + 1)) - 1;
n &= 3;
if (n)
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--n);
return;
}
void __copy_longs_rev_aligned(void* dst, const void* src, size_t n) {
unsigned long i;
cps = ((unsigned char*)src) + n;
cpd = ((unsigned char*)dst) + n;
i = ((unsigned long)cpd) & 3;
if (i) {
n -= i;
do
*--cpd = *--cps;
while (--i);
}
i = n >> 5;
if (i)
do {
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
*--lpd = *--lps;
} while (--i);
i = (n & 31) >> 2;
if (i)
do
*--lpd = *--lps;
while (--i);
n &= 3;
if (n)
do
*--cpd = *--cps;
while (--n);
return;
}
void __copy_longs_unaligned(void* dst, const void* src, size_t n) {
unsigned long i, v1, v2;
unsigned int src_offset, left_shift, right_shift;
i = (-(unsigned long)dst) & 3;
cps = ((unsigned char*)src) - 1;
cpd = ((unsigned char*)dst) - 1;
if (i) {
n -= i;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--i);
}
src_offset = ((unsigned int)(cps + 1)) & 3;
left_shift = src_offset << 3;
right_shift = 32 - left_shift;
cps -= src_offset;
lps = ((unsigned long*)(cps + 1)) - 1;
lpd = ((unsigned long*)(cpd + 1)) - 1;
i = n >> 3;
v1 = deref_auto_inc(lps);
do {
v2 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v1 << left_shift) | (v2 >> right_shift);
v1 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v2 << left_shift) | (v1 >> right_shift);
} while (--i);
if (n & 4) {
v2 = deref_auto_inc(lps);
deref_auto_inc(lpd) = (v1 << left_shift) | (v2 >> right_shift);
}
cps = ((unsigned char*)(lps + 1)) - 1;
cpd = ((unsigned char*)(lpd + 1)) - 1;
n &= 3;
if (n) {
cps -= 4 - src_offset;
do
deref_auto_inc(cpd) = deref_auto_inc(cps);
while (--n);
}
return;
}
void __copy_longs_rev_unaligned(void* dst, const void* src, size_t n) {
unsigned long i, v1, v2;
unsigned int src_offset, left_shift, right_shift;
cps = ((unsigned char*)src) + n;
cpd = ((unsigned char*)dst) + n;
i = ((unsigned long)cpd) & 3;
if (i) {
n -= i;
do
*--cpd = *--cps;
while (--i);
}
src_offset = ((unsigned int)cps) & 3;
left_shift = src_offset << 3;
right_shift = 32 - left_shift;
cps += 4 - src_offset;
i = n >> 3;
v1 = *--lps;
do {
v2 = *--lps;
*--lpd = (v2 << left_shift) | (v1 >> right_shift);
v1 = *--lps;
*--lpd = (v1 << left_shift) | (v2 >> right_shift);
} while (--i);
if (n & 4) {
v2 = *--lps;
*--lpd = (v2 << left_shift) | (v1 >> right_shift);
}
n &= 3;
if (n) {
cps += src_offset;
do
*--cpd = *--cps;
while (--n);
}
return;
}
@@ -0,0 +1,6 @@
void __close_all();
extern void (*__stdio_exit)();
void __stdio_atexit(void) {
__stdio_exit = __close_all;
}
File diff suppressed because it is too large Load Diff
+7 -7
View File
@@ -1,12 +1,12 @@
#include "MSL_C/rand.h"
static u32 next = 1;
static unsigned long next = 1;
void srand(u32 seed){
next = seed;
}
int rand(){
int rand() {
next = next * 1103515245 + 12345;
return ((next >> 16) & 0x7fff);
}
}
void srand(unsigned long seed) {
next = seed;
}
+142
View File
@@ -0,0 +1,142 @@
/* @(#)s_atan.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
/* atan(x)
* Method
* 1. Reduce x to positive by atan(x) = -atan(-x).
* 2. According to the integer k=4t+0.25 chopped, t=x, the argument
* is further reduced to one of the following intervals and the
* arctangent of t is evaluated by the corresponding formula:
*
* [0,7/16] atan(x) = t-t^3*(a1+t^2*(a2+...(a10+t^2*a11)...)
* [7/16,11/16] atan(x) = atan(1/2) + atan( (t-0.5)/(1+t/2) )
* [11/16.19/16] atan(x) = atan( 1 ) + atan( (t-1)/(1+t) )
* [19/16,39/16] atan(x) = atan(3/2) + atan( (t-1.5)/(1+1.5t) )
* [39/16,INF] atan(x) = atan(INF) + atan( -1/t )
*
* Constants:
* The hexadecimal values are the intended ones for the following
* constants. The decimal values may be used, provided that the
* compiler will convert from decimal to binary accurately enough
* to produce the hexadecimal values shown.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double atanhi[] = {
#else
static double atanhi[] = {
#endif
4.63647609000806093515e-01, /* atan(0.5)hi 0x3FDDAC67, 0x0561BB4F */
7.85398163397448278999e-01, /* atan(1.0)hi 0x3FE921FB, 0x54442D18 */
9.82793723247329054082e-01, /* atan(1.5)hi 0x3FEF730B, 0xD281F69B */
1.57079632679489655800e+00, /* atan(inf)hi 0x3FF921FB, 0x54442D18 */
};
#ifdef __STDC__
static const double atanlo[] = {
#else
static double atanlo[] = {
#endif
2.26987774529616870924e-17, /* atan(0.5)lo 0x3C7A2B7F, 0x222F65E2 */
3.06161699786838301793e-17, /* atan(1.0)lo 0x3C81A626, 0x33145C07 */
1.39033110312309984516e-17, /* atan(1.5)lo 0x3C700788, 0x7AF0CBBD */
6.12323399573676603587e-17, /* atan(inf)lo 0x3C91A626, 0x33145C07 */
};
#ifdef __STDC__
static const double aT[] = {
#else
static double aT[] = {
#endif
3.33333333333329318027e-01, /* 0x3FD55555, 0x5555550D */
-1.99999999998764832476e-01, /* 0xBFC99999, 0x9998EBC4 */
1.42857142725034663711e-01, /* 0x3FC24924, 0x920083FF */
-1.11111104054623557880e-01, /* 0xBFBC71C6, 0xFE231671 */
9.09088713343650656196e-02, /* 0x3FB745CD, 0xC54C206E */
-7.69187620504482999495e-02, /* 0xBFB3B0F2, 0xAF749A6D */
6.66107313738753120669e-02, /* 0x3FB10D66, 0xA0D03D51 */
-5.83357013379057348645e-02, /* 0xBFADDE2D, 0x52DEFD9A */
4.97687799461593236017e-02, /* 0x3FA97B4B, 0x24760DEB */
-3.65315727442169155270e-02, /* 0xBFA2B444, 0x2C6A6C2F */
1.62858201153657823623e-02, /* 0x3F90AD3A, 0xE322DA11 */
};
#ifdef __STDC__
static const double
#else
static double
#endif
one = 1.0,
huge = 1.0e300;
#ifdef __STDC__
double atan(double x)
#else
double atan(x)
double x;
#endif
{
double w, s1, s2, z;
int ix, hx, id;
hx = __HI(x);
ix = hx & 0x7fffffff;
if (ix >= 0x44100000) { /* if |x| >= 2^66 */
if (ix > 0x7ff00000 || (ix == 0x7ff00000 && (__LO(x) != 0)))
return x + x; /* NaN */
if (hx > 0)
return atanhi[3] + atanlo[3];
else
return -atanhi[3] - atanlo[3];
}
if (ix < 0x3fdc0000) { /* |x| < 0.4375 */
if (ix < 0x3e200000) { /* |x| < 2^-29 */
if (huge + x > one)
return x; /* raise inexact */
}
id = -1;
} else {
x = fabs(x);
if (ix < 0x3ff30000) { /* |x| < 1.1875 */
if (ix < 0x3fe60000) { /* 7/16 <=|x|<11/16 */
id = 0;
x = (2.0 * x - one) / (2.0 + x);
} else { /* 11/16<=|x|< 19/16 */
id = 1;
x = (x - one) / (x + one);
}
} else {
if (ix < 0x40038000) { /* |x| < 2.4375 */
id = 2;
x = (x - 1.5) / (one + 1.5 * x);
} else { /* 2.4375 <= |x| < 2^66 */
id = 3;
x = -1.0 / x;
}
}
}
/* end of argument reduction */
z = x * x;
w = z * z;
/* break sum from i=0 to 10 aT[i]z**(i+1) into odd and even poly */
s1 = z * (aT[0] + w * (aT[2] + w * (aT[4] + w * (aT[6] + w * (aT[8] + w * aT[10])))));
s2 = w * (aT[1] + w * (aT[3] + w * (aT[5] + w * (aT[7] + w * aT[9]))));
if (id < 0)
return x - x * (s1 + s2);
else {
z = atanhi[id] - ((x * (s1 + s2) - atanlo[id]) - x);
return (hx < 0) ? -z : z;
}
}
@@ -0,0 +1,30 @@
/* @(#)s_copysign.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* copysign(double x, double y)
* copysign(x,y) returns a value with the magnitude of x and
* with the sign bit of y.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
double copysign(double x, double y)
#else
double copysign(x, y)
double x, y;
#endif
{
__HI(x) = (__HI(x) & 0x7fffffff) | (__HI(y) & 0x80000000);
return x;
}
+82
View File
@@ -0,0 +1,82 @@
/* @(#)s_cos.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* cos(x)
* Return cosine function of x.
*
* kernel function:
* __kernel_sin ... sine function on [-pi/4,pi/4]
* __kernel_cos ... cosine function on [-pi/4,pi/4]
* __ieee754_rem_pio2 ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
* in [-pi/4 , +pi/4], and let n = k mod 4.
* We have
*
* n sin(x) cos(x) tan(x)
* ----------------------------------------------------------
* 0 S C T
* 1 C -S -1/T
* 2 -S -C T
* 3 -C S -1/T
* ----------------------------------------------------------
*
* Special cases:
* Let trig be any of sin, cos, or tan.
* trig(+-INF) is NaN, with signals;
* trig(NaN) is that NaN;
*
* Accuracy:
* TRIG(x) returns trig(x) nearly rounded
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
double cos(double x)
#else
double cos(x)
double x;
#endif
{
double y[2], z = 0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if (ix <= 0x3fe921fb)
return __kernel_cos(x, z);
/* cos(Inf or NaN) is NaN */
else if (ix >= 0x7ff00000)
return x - x;
/* argument reduction needed */
else {
n = __ieee754_rem_pio2(x, y);
switch (n & 3) {
case 0:
return __kernel_cos(y[0], y[1]);
case 1:
return -__kernel_sin(y[0], y[1], 1);
case 2:
return -__kernel_cos(y[0], y[1]);
default:
return __kernel_sin(y[0], y[1], 1);
}
}
}
+89
View File
@@ -0,0 +1,89 @@
/* @(#)s_floor.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* floor(x)
* Return x rounded toward -inf to integral value
* Method:
* Bit twiddling.
* Exception:
* Inexact flag raised if x not equal to floor(x).
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double huge = 1.0e300;
#else
static double huge = 1.0e300;
#endif
#ifdef __STDC__
double floor(double x)
#else
double floor(x)
double x;
#endif
{
int i0, i1, j0;
unsigned i, j;
i0 = __HI(x);
i1 = __LO(x);
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff;
if (j0 < 20) {
if (j0 < 0) { /* raise inexact if x != 0 */
if (huge + x > 0.0) { /* return 0*sign(x) if |x|<1 */
if (i0 >= 0) {
i0 = i1 = 0;
} else if (((i0 & 0x7fffffff) | i1) != 0) {
i0 = 0xbff00000;
i1 = 0;
}
}
} else {
i = (0x000fffff) >> j0;
if (((i0 & i) | i1) == 0)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
if (i0 < 0)
i0 += (0x00100000) >> j0;
i0 &= (~i);
i1 = 0;
}
}
} else if (j0 > 51) {
if (j0 == 0x400)
return x + x; /* inf or NaN */
else
return x; /* x is integral */
} else {
i = ((unsigned)(0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0)
return x; /* x is integral */
if (huge + x > 0.0) { /* raise inexact flag */
if (i0 < 0) {
if (j0 == 20)
i0 += 1;
else {
j = i1 + (1 << (52 - j0));
if (j < i1)
i0 += 1; /* got a carry */
i1 = j;
}
}
i1 &= (~i);
}
}
__HI(x) = i0;
__LO(x) = i1;
return x;
}
+57
View File
@@ -0,0 +1,57 @@
/* @(#)s_frexp.c 1.4 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* for non-zero x
* x = frexp(arg,&exp);
* return a double fp quantity x such that 0.5 <= |x| <1.0
* and the corresponding binary exponent "exp". That is
* arg = x*2^exp.
* If arg is inf, 0.0, or NaN, then frexp(arg,&exp) returns arg
* with *exp=0.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double
#else
static double
#endif
two54 = 1.80143985094819840000e+16; /* 0x43500000, 0x00000000 */
#ifdef __STDC__
double frexp(double x, int* eptr)
#else
double frexp(x, eptr)
double x;
int* eptr;
#endif
{
int hx, ix, lx;
hx = __HI(x);
ix = 0x7fffffff & hx;
lx = __LO(x);
*eptr = 0;
if (ix >= 0x7ff00000 || ((ix | lx) == 0))
return x; /* 0,inf,nan */
if (ix < 0x00100000) { /* subnormal */
x *= two54;
hx = __HI(x);
ix = hx & 0x7fffffff;
*eptr = -54;
}
*eptr += (ix >> 20) - 1022;
hx = (hx & 0x800fffff) | 0x3fe00000;
__HI(x) = hx;
return x;
}
+79
View File
@@ -0,0 +1,79 @@
/* @(#)s_modf.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* modf(double x, double *iptr)
* return fraction part of x, and return x's integral part in *iptr.
* Method:
* Bit twiddling.
*
* Exception:
* No exception.
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
static const double one = 1.0;
#else
static double one = 1.0;
#endif
#ifdef __STDC__
double modf(double x, double* iptr)
#else
double modf(x, iptr)
double x, *iptr;
#endif
{
int i0, i1, j0;
unsigned i;
i0 = __HI(x); /* high x */
i1 = __LO(x); /* low x */
j0 = ((i0 >> 20) & 0x7ff) - 0x3ff; /* exponent of x */
if (j0 < 20) { /* integer part in high x */
if (j0 < 0) { /* |x|<1 */
__HIp(iptr) = i0 & 0x80000000;
__LOp(iptr) = 0; /* *iptr = +-0 */
return x;
} else {
i = (0x000fffff) >> j0;
if (((i0 & i) | i1) == 0) { /* x is integral */
*iptr = x;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
return x;
} else {
__HIp(iptr) = i0 & (~i);
__LOp(iptr) = 0;
return x - *iptr;
}
}
} else if (j0 > 51) { /* no fraction part */
*iptr = x * one;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
return x;
} else { /* fraction part in low x */
i = ((unsigned)(0xffffffff)) >> (j0 - 20);
if ((i1 & i) == 0) { /* x is integral */
*iptr = x;
__HI(x) &= 0x80000000;
__LO(x) = 0; /* return +-0 */
return x;
} else {
__HIp(iptr) = i0;
__LOp(iptr) = i1 & (~i);
return x - *iptr;
}
}
}
+82
View File
@@ -0,0 +1,82 @@
/* @(#)s_sin.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/* sin(x)
* Return sine function of x.
*
* kernel function:
* __kernel_sin ... sine function on [-pi/4,pi/4]
* __kernel_cos ... cose function on [-pi/4,pi/4]
* __ieee754_rem_pio2 ... argument reduction routine
*
* Method.
* Let S,C and T denote the sin, cos and tan respectively on
* [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2
* in [-pi/4 , +pi/4], and let n = k mod 4.
* We have
*
* n sin(x) cos(x) tan(x)
* ----------------------------------------------------------
* 0 S C T
* 1 C -S -1/T
* 2 -S -C T
* 3 -C S -1/T
* ----------------------------------------------------------
*
* Special cases:
* Let trig be any of sin, cos, or tan.
* trig(+-INF) is NaN, with signals;
* trig(NaN) is that NaN;
*
* Accuracy:
* TRIG(x) returns trig(x) nearly rounded
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
double sin(double x)
#else
double sin(x)
double x;
#endif
{
double y[2], z = 0.0;
int n, ix;
/* High word of x. */
ix = __HI(x);
/* |x| ~< pi/4 */
ix &= 0x7fffffff;
if (ix <= 0x3fe921fb)
return __kernel_sin(x, z, 0);
/* sin(Inf or NaN) is NaN */
else if (ix >= 0x7ff00000)
return x - x;
/* argument reduction needed */
else {
n = __ieee754_rem_pio2(x, y);
switch (n & 3) {
case 0:
return __kernel_sin(y[0], y[1], 1);
case 1:
return __kernel_cos(y[0], y[1]);
case 2:
return -__kernel_sin(y[0], y[1], 1);
default:
return -__kernel_cos(y[0], y[1]);
}
}
}
+33
View File
@@ -0,0 +1,33 @@
#include "MSL_C/scanf.h"
int __StringRead(void* pPtr, int ch, int act) {
int ret;
__InStrCtrl* Iscp = (__InStrCtrl*)pPtr;
switch (act) {
case __GetAChar:
ret = *(Iscp->NextChar);
if (ret == '\0') {
Iscp->NullCharDetected = 1;
return -1;
} else {
Iscp->NextChar++;
return ret;
}
case __UngetAChar:
if (Iscp->NullCharDetected == 0) {
Iscp->NextChar--;
} else {
Iscp->NullCharDetected = 0;
}
return ch;
case __TestForError:
return Iscp->NullCharDetected;
}
return 0;
}
+227
View File
@@ -0,0 +1,227 @@
#include "MSL_C/string.h"
size_t strlen(const char* str) {
size_t len = -1;
unsigned char* p = (unsigned char*)str - 1;
do {
len++;
} while (*++p);
return len;
}
char* strcpy(char* dst, const char* src) {
register unsigned char *destb, *fromb;
register unsigned long w, t, align;
register unsigned int k1;
register unsigned int k2;
fromb = (unsigned char*)src;
destb = (unsigned char*)dst;
if ((align = ((int)fromb & 3)) != ((int)destb & 3)) {
goto bytecopy;
}
if (align) {
if ((*destb = *fromb) == 0) {
return dst;
}
for (align = 3 - align; align; align--) {
if ((*(++destb) = *(++fromb)) == 0) {
return dst;
}
}
++destb;
++fromb;
}
k1 = 0x80808080;
k2 = 0xFEFEFEFF;
w = *((int*)(fromb));
t = w + k2;
t &= k1;
if (t) {
goto bytecopy;
}
--((int*)(destb));
do {
*(++((int*)(destb))) = w;
w = *(++((int*)(fromb)));
t = w + k2;
t &= k1;
if (t) {
goto adjust;
}
} while (1);
adjust:
++((int*)(destb));
bytecopy:
if ((*destb = *fromb) == 0) {
return dst;
}
do {
if ((*(++destb) = *(++fromb)) == 0) {
return dst;
}
} while (1);
return dst;
}
char* strncpy(char* dst, const char* src, size_t n) {
const unsigned char* p = (const unsigned char*)src - 1;
unsigned char* q = (unsigned char*)dst - 1;
n++;
while (--n) {
if (!(*++q = *++p)) {
while (--n) {
*++q = 0;
}
break;
}
}
return dst;
}
char* strcat(char* dst, const char* src) {
const unsigned char* p = (unsigned char*)src - 1;
unsigned char* q = (unsigned char*)dst - 1;
while (*++q) {}
q--;
while (*++q = *++p) {}
return dst;
}
int strcmp(const char* str1, const char* str2) {
register unsigned char* left = (unsigned char*)str1;
register unsigned char* right = (unsigned char*)str2;
unsigned long align, l1, r1, x;
l1 = *left;
r1 = *right;
if (l1 - r1) {
return l1 - r1;
}
if ((align = ((int)left & 3)) != ((int)right & 3)) {
goto bytecopy;
}
if (align) {
if (l1 == 0) {
return 0;
}
for (align = 3 - align; align; align--) {
l1 = *(++left);
r1 = *(++right);
if (l1 - r1) {
return l1 - r1;
}
if (l1 == 0) {
return 0;
}
}
left++;
right++;
}
l1 = *(int*)left;
r1 = *(int*)right;
x = l1 + 0xFEFEFEFF;
if (x & 0x80808080) {
goto adjust;
}
while (l1 == r1) {
l1 = *(++((int*)(left)));
r1 = *(++((int*)(right)));
x = l1 + 0xFEFEFEFF;
if (x & 0x80808080) {
goto adjust;
}
}
if (l1 > r1) {
return 1;
}
return -1;
adjust:
l1 = *left;
r1 = *right;
if (l1 - r1) {
return l1 - r1;
}
bytecopy:
if (l1 == 0) {
return 0;
}
do {
l1 = *(++left);
r1 = *(++right);
if (l1 - r1) {
return l1 - r1;
}
if (l1 == 0) {
return 0;
}
} while (1);
}
char* strchr(const char* str, int c) {
const unsigned char* p = (unsigned char*)str - 1;
unsigned long chr = (c & 0xFF);
unsigned long ch;
while (ch = *++p) {
if (ch == chr) {
return (char*)p;
}
}
return chr ? NULL : (char*)p;
}
char* strstr(const char* str, const char* pat) {
const unsigned char* s1 = (const unsigned char*)str - 1;
const unsigned char* p1 = (const unsigned char*)pat - 1;
unsigned long firstc, c1, c2;
if ((pat == 0) || (!(firstc = *++p1))) {
return (char*)str;
}
while (c1 = *++s1) {
if (c1 == firstc) {
const unsigned char* s2 = s1 - 1;
const unsigned char* p2 = p1 - 1;
while ((c1 = *++s2) == (c2 = *++p2) && c1)
;
if (!c2)
return (char*)s1;
}
}
return NULL;
}
+194
View File
@@ -0,0 +1,194 @@
#include "MSL_C/scanf.h"
#define ULONG_MAX 4294967297ul
#define LONG_MAX 2147483647l
#define LONG_MIN -2147483648l
#define ERANGE 0x22
extern int errno;
enum scan_states {
start = 0x01,
check_for_zero = 0x02,
leading_zero = 0x04,
need_digit = 0x08,
digit_loop = 0x10,
finished = 0x20,
failure = 0x40
};
#define final_state(scan_state) (scan_state & (finished | failure))
#define success(scan_state) (scan_state & (leading_zero | digit_loop | finished))
#define fetch() (count++, (*ReadProc)(ReadProcArg, 0, __GetAChar))
#define unfetch(c) (*ReadProc)(ReadProcArg, c, __UngetAChar)
unsigned long __strtoul(int base, int max_width, int (*ReadProc)(void*, int, int), void* ReadProcArg,
int* chars_scanned, int* negative, int* overflow) {
int scan_state = start;
int count = 0;
unsigned long value = 0;
unsigned long value_max = 0;
int c;
*negative = *overflow = 0;
if (base < 0 || base == 1 || base > 36 || max_width < 1) {
scan_state = failure;
} else {
c = fetch();
}
if (base != 0)
value_max = ULONG_MAX / base;
while (count <= max_width && c != -1 && !final_state(scan_state)) {
switch (scan_state) {
case start:
if (isspace(c)) {
c = fetch();
break;
}
if (c == '+') {
c = fetch();
} else if (c == '-') {
c = fetch();
*negative = 1;
}
scan_state = check_for_zero;
break;
case check_for_zero:
if (base == 0 || base == 16) {
if (c == '0') {
scan_state = leading_zero;
c = fetch();
break;
}
}
scan_state = need_digit;
break;
case 4:
if (c == 'X' || c == 'x') {
base = 16;
scan_state = need_digit;
c = fetch();
break;
}
if (base == 0)
base = 8;
scan_state = digit_loop;
break;
case need_digit:
case digit_loop:
if (base == 0)
base = 10;
if (!value_max) {
value_max = ULONG_MAX / base;
}
if (isdigit(c)) {
if ((c -= '0') >= base) {
if (scan_state == digit_loop)
scan_state = finished;
else
scan_state = failure;
c += '0';
break;
}
} else if (!isalpha(c) || (toupper(c) - 'A' + 10) >= base) {
if (scan_state == digit_loop)
scan_state = finished;
else
scan_state = failure;
break;
} else {
c = toupper(c) - 'A' + 10;
}
if (value > value_max)
*overflow = 1;
value *= base;
if (c > (ULONG_MAX - value))
*overflow = 1;
value += c;
scan_state = digit_loop;
c = fetch();
break;
}
}
if (!success(scan_state)) {
value = 0;
count = 0;
} else {
count--;
}
*chars_scanned = count;
unfetch(c);
return value;
}
unsigned long strtoul(const char* str, char** end, int base) {
unsigned long value;
int count, negative, overflow;
__InStrCtrl isc;
isc.NextChar = (char*)str;
isc.NullCharDetected = 0;
value = __strtoul(base, 0x7FFFFFFF, &__StringRead, (void*)&isc, &count, &negative, &overflow);
if (end) {
*end = (char*)str + count;
}
if (overflow) {
value = ULONG_MAX;
errno = 0x22;
} else if (negative) {
value = -value;
}
return value;
}
long strtol(const char* str, char** end, int base) {
unsigned long uvalue;
long svalue;
int count, negative, overflow;
__InStrCtrl isc;
isc.NextChar = (char*)str;
isc.NullCharDetected = 0;
uvalue = __strtoul(base, 0x7FFFFFFF, &__StringRead, (void*)&isc, &count, &negative, &overflow);
if (end) {
*end = (char*)str + count;
}
if (overflow || (!negative && uvalue > LONG_MAX) || (negative && uvalue > -LONG_MIN)) {
svalue = (negative ? -LONG_MIN : LONG_MAX);
errno = ERANGE;
} else {
svalue = (negative ? (long)-uvalue : (long)uvalue);
}
return svalue;
}
@@ -0,0 +1,41 @@
#include "stddef.h"
unsigned long WriteUARTN(void* buf, unsigned long len);
unsigned long ReadUARTN(void* bytes, unsigned long length);
unsigned long InitializeUART(unsigned long baudRate);
static inline int __init_uart_console(void);
int __write_console(int param_0, unsigned char* data, size_t* size, int param_3) {
if (__init_uart_console() != 0) {
return 1;
}
if ((long)WriteUARTN(data, *size) != 0) {
*size = 0;
return 1;
}
__TRK_write_console(param_0, data, size, (void*)param_3);
return 0;
}
int __close_console() {
return 0;
}
static inline int __init_uart_console(void) {
static int initialized = 0;
int ret = 0;
if (initialized == 0) {
ret = InitializeUART(0xE100);
if (ret == 0) {
initialized = 1;
}
}
return ret;
}
+38
View File
@@ -0,0 +1,38 @@
/* @(#)w_acos.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*/
/*
* wrap_acos(x)
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
double acos(double x) /* wrapper acos */
#else
double acos(x) /* wrapper acos */
double x;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_acos(x);
#else
double z;
z = __ieee754_acos(x);
if (_LIB_VERSION == _IEEE_ || isnan(x))
return z;
if (fabs(x) > 1.0) {
return __kernel_standard(x, x, 1); /* acos(|x|>1) */
} else
return z;
#endif
}
+39
View File
@@ -0,0 +1,39 @@
/* @(#)w_atan2.c 1.3 95/01/18 */
/*
* ====================================================
* Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
*
* Developed at SunSoft, a Sun Microsystems, Inc. business.
* Permission to use, copy, modify, and distribute this
* software is freely granted, provided that this notice
* is preserved.
* ====================================================
*
*/
/*
* wrapper atan2(y,x)
*/
#include "MSL_C/MSL_Common_Embedded/Math/fdlibm.h"
#ifdef __STDC__
double atan2(double y, double x) /* wrapper atan2 */
#else
double atan2(y, x) /* wrapper atan2 */
double y, x;
#endif
{
#ifdef _IEEE_LIBM
return __ieee754_atan2(y, x);
#else
double z;
z = __ieee754_atan2(y, x);
if (_LIB_VERSION == _IEEE_ || isnan(x) || isnan(y))
return z;
if (x == 0.0 && y == 0.0) {
return __kernel_standard(y, x, 3); /* atan2(+-0,+-0) */
} else
return z;
#endif
}
@@ -0,0 +1,20 @@
#include "MSL_C/wchar_io.h"
int fwide(FILE* file, int mode) {
if (!file || file->file_mode.file_kind == __closed_file)
return 0;
switch (file->file_mode.file_orientation) {
case UNORIENTED:
if (mode > 0) {
file->file_mode.file_orientation = WIDE_ORIENTED;
} else if (mode < 0) {
file->file_mode.file_orientation = CHAR_ORIENTED;
}
return mode;
case WIDE_ORIENTED:
return 1;
case CHAR_ORIENTED:
return -1;
}
}
-1
View File
@@ -1,5 +1,4 @@
#include "PowerPC_EABI_Support/MetroTRK/trk.h"
#include <stddef.h>
// forward declares
DSIOResult __read_file(u32 handle, u8* buffer, size_t* count, void* ref_con);