Home | History | Annotate | Line # | Download | only in fpsp 
      1 *	$NetBSD: scosh.sa,v 1.2 1994/10/26 07:49:39 cgd Exp $
      2 
      3 *	MOTOROLA MICROPROCESSOR & MEMORY TECHNOLOGY GROUP
      4 *	M68000 Hi-Performance Microprocessor Division
      5 *	M68040 Software Package 
      6 *
      7 *	M68040 Software Package Copyright (c) 1993, 1994 Motorola Inc.
      8 *	All rights reserved.
      9 *
     10 *	THE SOFTWARE is provided on an "AS IS" basis and without warranty.
     11 *	To the maximum extent permitted by applicable law,
     12 *	MOTOROLA DISCLAIMS ALL WARRANTIES WHETHER EXPRESS OR IMPLIED,
     13 *	INCLUDING IMPLIED WARRANTIES OF MERCHANTABILITY OR FITNESS FOR A
     14 *	PARTICULAR PURPOSE and any warranty against infringement with
     15 *	regard to the SOFTWARE (INCLUDING ANY MODIFIED VERSIONS THEREOF)
     16 *	and any accompanying written materials. 
     17 *
     18 *	To the maximum extent permitted by applicable law,
     19 *	IN NO EVENT SHALL MOTOROLA BE LIABLE FOR ANY DAMAGES WHATSOEVER
     20 *	(INCLUDING WITHOUT LIMITATION, DAMAGES FOR LOSS OF BUSINESS
     21 *	PROFITS, BUSINESS INTERRUPTION, LOSS OF BUSINESS INFORMATION, OR
     22 *	OTHER PECUNIARY LOSS) ARISING OF THE USE OR INABILITY TO USE THE
     23 *	SOFTWARE.  Motorola assumes no responsibility for the maintenance
     24 *	and support of the SOFTWARE.  
     25 *
     26 *	You are hereby granted a copyright license to use, modify, and
     27 *	distribute the SOFTWARE so long as this entire notice is retained
     28 *	without alteration in any modified and/or redistributed versions,
     29 *	and that such modified versions are clearly identified as such.
     30 *	No licenses are granted by implication, estoppel or otherwise
     31 *	under any patents or trademarks of Motorola, Inc.
     32 
     33 *
     34 *	scosh.sa 3.1 12/10/90
     35 *
     36 *	The entry point sCosh computes the hyperbolic cosine of
     37 *	an input argument; sCoshd does the same except for denormalized
     38 *	input.
     39 *
     40 *	Input: Double-extended number X in location pointed to
     41 *		by address register a0.
     42 *
     43 *	Output: The value cosh(X) returned in floating-point register Fp0.
     44 *
     45 *	Accuracy and Monotonicity: The returned result is within 3 ulps in
     46 *		64 significant bit, i.e. within 0.5001 ulp to 53 bits if the
     47 *		result is subsequently rounded to double precision. The
     48 *		result is provably monotonic in double precision.
     49 *
     50 *	Speed: The program sCOSH takes approximately 250 cycles.
     51 *
     52 *	Algorithm:
     53 *
     54 *	COSH
     55 *	1. If |X| > 16380 log2, go to 3.
     56 *
     57 *	2. (|X| <= 16380 log2) Cosh(X) is obtained by the formulae
     58 *		y = |X|, z = exp(Y), and
     59 *		cosh(X) = (1/2)*( z + 1/z ).
     60 *		Exit.
     61 *
     62 *	3. (|X| > 16380 log2). If |X| > 16480 log2, go to 5.
     63 *
     64 *	4. (16380 log2 < |X| <= 16480 log2)
     65 *		cosh(X) = sign(X) * exp(|X|)/2.
     66 *		However, invoking exp(|X|) may cause premature overflow.
     67 *		Thus, we calculate sinh(X) as follows:
     68 *		Y	:= |X|
     69 *		Fact	:=	2**(16380)
     70 *		Y'	:= Y - 16381 log2
     71 *		cosh(X) := Fact * exp(Y').
     72 *		Exit.
     73 *
     74 *	5. (|X| > 16480 log2) sinh(X) must overflow. Return
     75 *		Huge*Huge to generate overflow and an infinity with
     76 *		the appropriate sign. Huge is the largest finite number in
     77 *		extended format. Exit.
     78 *
     79 
     80 SCOSH	IDNT	2,1 Motorola 040 Floating Point Software Package
     81 
     82 	section	8
     83 
     84 	xref	t_ovfl
     85 	xref	t_frcinx
     86 	xref	setox
     87 
     88 T1	DC.L $40C62D38,$D3D64634 ... 16381 LOG2 LEAD
     89 T2	DC.L $3D6F90AE,$B1E75CC7 ... 16381 LOG2 TRAIL
     90 
     91 TWO16380 DC.L $7FFB0000,$80000000,$00000000,$00000000
     92 
     93 	xdef	scoshd
     94 scoshd:
     95 *--COSH(X) = 1 FOR DENORMALIZED X
     96 
     97 	FMOVE.S		#:3F800000,FP0
     98 
     99 	FMOVE.L		d1,FPCR
    100 	FADD.S		#:00800000,FP0
    101 	bra		t_frcinx
    102 
    103 	xdef	scosh
    104 scosh:
    105 	FMOVE.X		(a0),FP0	...LOAD INPUT
    106 
    107 	move.l		(a0),d0
    108 	move.w		4(a0),d0
    109 	ANDI.L		#$7FFFFFFF,d0
    110 	CMPI.L		#$400CB167,d0
    111 	BGT.B		COSHBIG
    112 
    113 *--THIS IS THE USUAL CASE, |X| < 16380 LOG2
    114 *--COSH(X) = (1/2) * ( EXP(X) + 1/EXP(X) )
    115 
    116 	FABS.X		FP0		...|X|
    117 
    118 	move.l		d1,-(sp)
    119 	clr.l		d1
    120 	fmovem.x	fp0,(a0)	;pass parameter to setox
    121 	bsr		setox		...FP0 IS EXP(|X|)
    122 	FMUL.S		#:3F000000,FP0	...(1/2)EXP(|X|)
    123 	move.l		(sp)+,d1
    124 
    125 	FMOVE.S		#:3E800000,FP1	...(1/4)
    126 	FDIV.X		FP0,FP1	 	...1/(2 EXP(|X|))
    127 
    128 	FMOVE.L		d1,FPCR
    129 	FADD.X		fp1,FP0
    130 
    131 	bra		t_frcinx
    132 
    133 COSHBIG:
    134 	CMPI.L		#$400CB2B3,d0
    135 	BGT.B		COSHHUGE
    136 
    137 	FABS.X		FP0
    138 	FSUB.D		T1(pc),FP0		...(|X|-16381LOG2_LEAD)
    139 	FSUB.D		T2(pc),FP0		...|X| - 16381 LOG2, ACCURATE
    140 
    141 	move.l		d1,-(sp)
    142 	clr.l		d1
    143 	fmovem.x	fp0,(a0)
    144 	bsr		setox
    145 	fmove.l		(sp)+,fpcr
    146 
    147 	FMUL.X		TWO16380(pc),FP0
    148 	bra		t_frcinx
    149 
    150 COSHHUGE:
    151 	fmove.l		#0,fpsr		;clr N bit if set by source
    152 	bclr.b		#7,(a0)		;always return positive value
    153 	fmovem.x	(a0),fp0
    154 	bra		t_ovfl
    155 
    156 	end
    157