diff --git a/fenv.asm b/fenv.asm deleted file mode 100644 index 843f063..0000000 --- a/fenv.asm +++ /dev/null @@ -1,361 +0,0 @@ - keep obj/fenv - mcopy fenv.macros - case on - -**************************************************************** -* -* Fenv - Floating-point environment access -* -* This code provides routines to query and modify the -* floating-point environment. -* -* Note: This relies on and only works with SANE. -* -**************************************************************** -* -fenv private dummy segment - end - -FE_ALL_EXCEPT gequ $001F - -**************************************************************** -* -* int feclearexcept(int excepts); -* -* Clear floating-point exceptions -* -* Inputs: -* excepts - floating-point exceptions to clear -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -feclearexcept start - - csubroutine (2:excepts),0 - - FGETENV get current environment - phx - - lda excepts - and #FE_ALL_EXCEPT - eor #$FFFF mask off excepts to clear - xba - and 1,S - sta 1,S - FSETENV clear them - - creturn 2:#0 - end - -**************************************************************** -* -* int fegetexceptflag(fexcept_t *flagp, int excepts); -* -* Get floating-point exception flags. -* -* Inputs: -* flagp - pointer to location to store exception flags -* excepts - floating-point exceptions to get -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -fegetexceptflag start - - csubroutine (4:flagp,2:excepts),0 - - FGETENV get current environment - tya - and excepts get desired exceptions - and #FE_ALL_EXCEPT - sta [flagp] store them in *flagp - - creturn 2:#0 - end - -**************************************************************** -* -* int feraiseexcept(int excepts); -* -* Raise floating-point exceptions -* -* Inputs: -* excepts - floating-point exceptions to raise -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -feraiseexcept start - - csubroutine (2:excepts),0 - - lda excepts - and #FE_ALL_EXCEPT - beq done - pha - FSETXCP raise exceptions - -done creturn 2:#0 - end - -**************************************************************** -* -* int fesetexceptflag(fexcept_t *flagp, int excepts); -* -* Set (but do not raise) floating-point exception flags -* -* Inputs: -* flagp - pointer to stored exception flags -* excepts - floating-point exceptions to set -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** - -fesetexceptflag start - - csubroutine (4:flagp,2:excepts),0 - - FGETENV get env with excepts masked off - phx - lda excepts - and #FE_ALL_EXCEPT - eor #$FFFF - xba - and 1,S - sta 1,S - - lda [flagp] set new exceptions - and excepts - and #FE_ALL_EXCEPT - xba - ora 1,S - sta 1,S - FSETENV - - creturn 2:#0 - end - -**************************************************************** -* -* int fetestexcept(int excepts); -* -* Test if floating-point exception flags are set -* -* Inputs: -* excepts - floating-point exceptions to test for -* -* Outputs: -* Bitwise or of exceptions that are set -* -**************************************************************** -* -fetestexcept start - - csubroutine (2:excepts),0 - - FGETENV get exception flags - tya - and excepts mask to just the ones we want - and #FE_ALL_EXCEPT - sta excepts - - creturn 2:excepts - end - -**************************************************************** -* -* int fegetround(void); -* -* Get the current rounding direction -* -* Outputs: -* The current rounding direction -* -**************************************************************** -* -fegetround start - FGETENV get high word of environment - tya - and #$00C0 just rounding direction - rtl - end - -**************************************************************** -* -* int fesetround(int round); -* -* Set the current rounding direction -* -* Inputs: -* round - the rounding direction to set -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -fesetround start - - csubroutine (2:round),0 - - lda round flip words - xba - sta round - and #$3FFF do nothing if not a valid rounding dir - bne done - - FGETENV set the rounding direction - txa - and #$3FFF - ora round - pha - FSETENV - - stz round -done creturn 2:round - end - -**************************************************************** -* -* int fegetenv(fenv_t *envp); -* -* Get the current floating-point environment -* -* Inputs: -* envp - pointer to location to store environment -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -fegetenv start - - csubroutine (4:envp),0 - - FGETENV get the environment - txa - sta [envp] store it in *envp - - creturn 2:#0 - end - -**************************************************************** -* -* int feholdexcept(fenv_t *envp); -* -* Get environment, then clear status flags and disable halts -* -* Inputs: -* envp - pointer to location to store environment -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -feholdexcept start - - csubroutine (4:envp),0 - - FGETENV get the environment - txa - sta [envp] store it in *envp - - and #$E0E0 clear exception flags and disable halts - pha - FSETENV set the new environment - - creturn 2:#0 - end - -**************************************************************** -* -* int fesetenv(const fenv_t *envp); -* -* Set the floating-point environment -* -* Inputs: -* envp - pointer to environment to set -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -fesetenv start - - csubroutine (4:envp),0 - - lda [envp] set the environment - pha - FSETENV - - creturn 2:#0 - end - -**************************************************************** -* -* int feupdateenv(const fenv_t *envp); -* -* Save exceptions, set environment, then re-raise exceptions -* -* Inputs: -* envp - pointer to environment to set -* -* Outputs: -* Returns 0 if successful, non-zero otherwise -* -**************************************************************** -* -feupdateenv start - - csubroutine (4:envp),0 - - lda [envp] set the environment - pha - FPROCEXIT - - creturn 2:#0 - end - -**************************************************************** -* -* Default floating-point environment -* -**************************************************************** -* -__FE_DFL_ENV start - dc i2'0' - end - -**************************************************************** -* -* int __get_flt_rounds(void); -* -* Get the value of FLT_ROUNDS, accounting for rounding mode -* -* Outputs: -* Current value of FLT_ROUNDS -* -**************************************************************** -* -__get_flt_rounds start - FGETENV - tya get rounding direction in low bits of A - asl a - asl a - xba - inc a convert to values used by FLT_ROUNDS - and #$0003 - rtl - end diff --git a/fenv.macros b/fenv.macros deleted file mode 100644 index aa38744..0000000 --- a/fenv.macros +++ /dev/null @@ -1,117 +0,0 @@ - MACRO -&lab csubroutine &parms,&work -&lab anop - aif c:&work,.a - lclc &work -&work setc 0 -.a - gbla &totallen - gbla &worklen -&worklen seta &work -&totallen seta 0 - aif c:&parms=0,.e - lclc &len - lclc &p - lcla &i -&i seta 1 -.b -&p setc &parms(&i) -&len amid &p,2,1 - aif "&len"=":",.c -&len amid &p,1,2 -&p amid &p,4,l:&p-3 - ago .d -.c -&len amid &p,1,1 -&p amid &p,3,l:&p-2 -.d -&p equ &totallen+4+&work -&totallen seta &totallen+&len -&i seta &i+1 - aif &i<=c:&parms,^b -.e - tsc - aif &work=0,.f - sec - sbc #&work - tcs -.f - phd - tcd - mend - MACRO -&lab creturn &r -&lab anop - lclc &len - aif c:&r,.a - lclc &r -&r setc 0 -&len setc 0 - ago .h -.a -&len amid &r,2,1 - aif "&len"=":",.b -&len amid &r,1,2 -&r amid &r,4,l:&r-3 - ago .c -.b -&len amid &r,1,1 -&r amid &r,3,l:&r-2 -.c - aif &len<>2,.d - ldy &r - ago .h -.d - aif &len<>4,.e - ldx &r+2 - ldy &r - ago .h -.e - aif &len<>10,.g - ldy #&r - ldx #^&r - ago .h -.g - mnote 'Not a valid return length',16 - mexit -.h - aif &totallen=0,.i - lda &worklen+2 - sta &worklen+&totallen+2 - lda &worklen+1 - sta &worklen+&totallen+1 -.i - pld - tsc - clc - adc #&worklen+&totallen - tcs - aif &len=0,.j - tya -.j - rtl - mend - MACRO -&LAB FGETENV -&LAB PEA $03 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSETENV -&LAB PEA $01 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSETXCP -&LAB PEA $15 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FPROCEXIT -&LAB PEA $19 - LDX #$090A - JSL $E10000 - MEND diff --git a/fpextra.asm b/fpextra.asm deleted file mode 100644 index 0c24d5e..0000000 --- a/fpextra.asm +++ /dev/null @@ -1,116 +0,0 @@ - keep obj/fpextra - mcopy fpextra.macros - -**************************************************************** -* -* FPextra - extra floating-point routines -* -* This code provides routines dealing with floating-point -* numbers that are used only by ORCA/C, supplementing the -* ones in SysFloat. -* -**************************************************************** -* -fpextra private dummy segment - end - -**************************************************************** -* -* ~SinglePrecision - limit fp value to single precision & range -* -* Inputs: -* extended-format real on stack -* -**************************************************************** -* -~SinglePrecision start - tsc - clc - adc #4 - ldy #0 - phy - pha - phy - pha - phy - pha - phy - pha - FX2S - FS2X - rtl - end - -**************************************************************** -* -* ~DoublePrecision - limit fp value to double precision & range -* -* Inputs: -* extended-format real on stack -* -**************************************************************** -* -~DoublePrecision start - tsc - clc - adc #4 - ldy #0 - phy - pha - phy - pha - phy - pha - phy - pha - FX2D - FD2X - rtl - end - -**************************************************************** -* -* ~CompPrecision - limit fp value to comp precision & range -* -* Inputs: -* extended-format real on stack -* -* Note: This avoids calling FX2C on negative numbers, -* because it is buggy for certain values. -* -**************************************************************** -* -~CompPrecision start - tsc round to integer - clc - adc #4 - pea 0 - pha - FRINTX - lda 4+8,s - pha save original sign - asl a force sign to positive - lsr a - sta 6+8,s - tsc limit precision - clc - adc #6 - ldy #0 - phy - pha - phy - pha - phy - pha - phy - pha - FX2C - FC2X - pla restore original sign - bpl ret - lda 4+8,s - ora #$8000 - sta 4+8,s -ret rtl - end - diff --git a/fpextra.macros b/fpextra.macros deleted file mode 100644 index bda9964..0000000 --- a/fpextra.macros +++ /dev/null @@ -1,42 +0,0 @@ - MACRO -&LAB FX2S -&LAB PEA $0210 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2D -&LAB PEA $0110 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2C -&LAB PEA $0510 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FC2X -&LAB PEA $050E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FD2X -&LAB PEA $010E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FS2X -&LAB PEA $020E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FRINTX -&LAB PEA $0014 - LDX #$090A - JSL $E10000 - MEND diff --git a/int64.asm b/int64.asm index 297c85b..f5c5961 100644 --- a/int64.asm +++ b/int64.asm @@ -515,273 +515,3 @@ loop1 asl num1 do the remaining shift rt0 pld rtl rtl end - -**************************************************************** -* -* ~CnvULongLongReal - convert an unsigned long long integer -* into an extended SANE real -* -* Inputs: -* unsigned long long int on stack -* -* Outputs: -* extended real on stack -* -**************************************************************** -* -~CnvULongLongReal start -mantissa equ 4 mantissa (integer and fraction) -exponent equ mantissa+8 biased exponent and sign bit - - lda 1,S move return value - pha - lda 4,S - sta 2,S - tsc set up DP - phd - tcd - - lda mantissa+2 move 64-bit value to mantissa - sta mantissa - lda mantissa+4 - sta mantissa+2 - lda mantissa+6 - sta mantissa+4 - lda mantissa+8 - sta mantissa+6 - - ora mantissa if value is 0 then - ora mantissa+2 - ora mantissa+4 - beq ret return - - lda #63+16383 set initial exponent (2^63) and sign - sta exponent - - lda mantissa+6 if number is normalized (i=1) then - bmi ret return - -lp1 dec exponent normalize number - asl mantissa - rol mantissa+2 - rol mantissa+4 - rol mantissa+6 - bpl lp1 - -ret pld - rtl - end - -**************************************************************** -* -* ~CnvLongLongReal - convert a long long integer into -* an extended SANE real -* -* Inputs: -* signed long long int on stack -* -* Outputs: -* extended real on stack -* -**************************************************************** -* -~CnvLongLongReal start -mantissa equ 4 mantissa (integer and fraction) -exponent equ mantissa+8 biased exponent and sign bit - - lda 1,S move return value - pha - lda 4,S - sta 2,S - tsc set up DP - phd - tcd - - lda mantissa+2 move 64-bit value to mantissa - sta mantissa - lda mantissa+4 - sta mantissa+2 - lda mantissa+6 - sta mantissa+4 - lda mantissa+8 - sta mantissa+6 - - ora mantissa if value is 0 then - ora mantissa+2 - ora mantissa+4 - beq ret return - - ldy #0 default sign bit is 0 (positive) - lda mantissa+6 if mantissa is negative then - bpl lb0 - negate8 mantissa negate it - ldy #$8000 sign bit is 1 (negative) - -lb0 tya set sign - ora #63+16383 set initial exponent (2^63) - sta exponent - - lda mantissa+6 if number is normalized (i=1) then - bmi ret return - -lp1 dec exponent normalize number - asl mantissa - rol mantissa+2 - rol mantissa+4 - rol mantissa+6 - bpl lp1 - -ret pld - rtl - end - -**************************************************************** -* -* ~CnvRealLongLong - convert an extended SANE real into -* a long long integer -* -* Inputs: -* extended real on stack -* -* Outputs: -* signed long long int on stack -* -* Note: This avoids calling FX2C on negative numbers, -* because it is buggy for certain values. -* -**************************************************************** -* -~CnvRealLongLong start - tsc - clc - adc #4 - pea 0 push src address for fcpxx - pha - pea llmin|-16 push dst address for fcpxx - pea llmin - pea 0 push operand address for ftintx - pha - ftintx round - fcpxx compare with LLONG_MIN - bne convert - - lda #$8000 if it is LONG_MIN, use that value - sta 12,s - asl a - sta 10,s - sta 8,s - sta 6,s - bra done otherwise - -convert lda 4+8,s - pha save original sign - asl a force sign to positive - lsr a - sta 6+8,s - tsc - clc - adc #6 - pea 0 push src address for fx2c - pha - pea 0 push dst address for fx2c - inc a - inc a - pha - fx2c convert - pla if original value was negative - bpl done - sec - ldx #0 negate result - txa - sbc 6,s - sta 6,s - txa - sbc 6+2,s - sta 6+2,s - txa - sbc 6+4,s - sta 6+4,s - txa - sbc 6+6,s - sta 6+6,s - -done phb move return address - pla - plx - ply - phx - pha - plb - rtl - -llmin dc e'-9223372036854775808' - end - -**************************************************************** -* -* ~CnvRealULongLong - convert an extended SANE real into -* an unsigned long long integer -* -* Inputs: -* extended real on stack -* -* Outputs: -* unsigned long long int on stack -* -**************************************************************** -* -~CnvRealULongLong start - pea 0 initially assume val <= LLONG_MAX - - tsc - clc - adc #6 - pea 0 push src address for fcpxx - pha - pea llbig|-16 push dst address for fcpxx - pea llbig - pea 0 push operand address for ftintx - pha - ftintx round - fcpxx compare with LLONG_MAX+1 - bmi convert - - lda #1 if val > LLONG_MAX: - sta 1,S save flag to indicate this - tsc - clc - adc #6 - pea llbig|-16 push src address for fsubx - pea llbig - pea 0 push dst address for fsubx - pha - fsubx val -= LLONG_MAX+1 - -convert tsc - clc - adc #6 - pea 0 push src address for fx2c - pha - pea 0 push dst address for fx2c - inc a - inc a - pha - fx2c convert val as comp - - pla if orig val was > LLONG_MAX: - beq done - lda 12,s - eor #$8000 - sta 12,s result += LLONG_MAX+1 - -done phb move return address - pla - plx - ply - phx - pha - plb - rtl - -llbig dc e'9223372036854775808' - end diff --git a/int64.macros b/int64.macros index f4da0ae..aedead3 100644 --- a/int64.macros +++ b/int64.macros @@ -1,23 +1,4 @@ macro -&l negate8 &n1 -&l ~setm - sec - ldy #0 - tya - sbc &n1 - sta &n1 - tya - sbc &n1+2 - sta &n1+2 - tya - sbc &n1+4 - sta &n1+4 - tya - sbc &n1+6 - sta &n1+6 - ~restm - mend - macro &l move4 &m1,&m2 lclb &yistwo &l ~setm @@ -140,27 +121,3 @@ .d sta 2+&op mend - MACRO -&LAB FTINTX -&LAB PEA $0016 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2C -&LAB PEA $0510 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCPXX -&LAB PEA $0A - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSUBX -&LAB PEA 2 - LDX #$090A - JSL $E10000 - MEND diff --git a/make b/make index 0976585..3e2ef86 100644 --- a/make +++ b/make @@ -19,7 +19,7 @@ if {#} == 0 unset exit end - for i in cc ctype string stdlib time setjmp orca fcntl vars toolglue signal int64 fenv fpextra math2 locale uchar + for i in cc ctype string stdlib time setjmp orca fcntl vars toolglue signal int64 locale uchar Newer obj/{i}.a {i}.asm if {Status} != 0 set exit on @@ -40,7 +40,7 @@ delete orcalib set list vars.a assert.a cc.a setjmp.a ctype.a string.a stdlib.a set list {list} time.a signal.a toolglue.a orca.a fcntl.a stdio.a int64.a -set list {list} fenv.a fpextra.a math2.a locale.a uchar.a +set list {list} locale.a uchar.a for i in {list} echo makelib orcalib +obj/{i} makelib orcalib +obj/{i} diff --git a/math2.asm b/math2.asm deleted file mode 100644 index 1cfc1f9..0000000 --- a/math2.asm +++ /dev/null @@ -1,4257 +0,0 @@ - keep obj/math2 - mcopy math2.macros - case on - -**************************************************************** -* -* Math2 - additional math routines -* -* This code provides additional functions from -* (including internal helper functions used by macros), -* supplementing the ones in SysFloat. -* -**************************************************************** - -math2 private dummy segment - copy equates.asm - end - -INVALID gequ $0001 exceptions -UNDERFLOW gequ $0002 -OVERFLOW gequ $0004 -DIVBYZERO gequ $0008 -INEXACT gequ $0010 - -TONEAREST gequ 0 rounding directions -UPWARD gequ 1 -DOWNWARD gequ 2 -TOWARDZERO gequ 3 - -**************************************************************** -* -* MathCommon2 - common work areas for the math library -* -**************************************************************** -* -MathCommon2 privdata -; -; temporary work space/return value -; -t1 ds 10 - end - -**************************************************************** -* -* int __fpclassifyf(float x); -* -* Classify a float value -* -* Inputs: -* val - the number to classify -* -* Outputs: -* one of the FP_* classification values -* -**************************************************************** -* -__fpclassifyf start - - csubroutine (10:val),0 - - tdc - clc - adc #val - ldy #0 - phy - pha - phy - pha - phy - pha - FX2S - FCLASSS - txa - and #$00FF - cmp #$00FC - bne lb1 - inc a -lb1 sta val - - creturn 2:val - end - -**************************************************************** -* -* int __fpclassifyd(double x); -* -* Classify a double value -* -* Inputs: -* val - the number to classify -* -* Outputs: -* one of the FP_* classification values -* -**************************************************************** -* -__fpclassifyd start - - csubroutine (10:val),0 - - tdc - clc - adc #val - ldy #0 - phy - pha - phy - pha - phy - pha - FX2D - FCLASSD - txa - and #$00FF - cmp #$00FC - bne lb1 - inc a -lb1 sta val - - creturn 2:val - end - -**************************************************************** -* -* int __fpclassifyl(long double x); -* -* Classify a long double value -* -* Inputs: -* val - the number to classify -* -* Outputs: -* one of the FP_* classification values -* -**************************************************************** -* -__fpclassifyl start - - csubroutine (10:val),0 - - tdc - clc - adc #val - pea 0 - pha - FCLASSX - txa - and #$00FF - cmp #$00FC - bne lb1 - inc a -lb1 sta val - - creturn 2:val - end - -**************************************************************** -* -* int __signbit(long double x); -* -* Get the sign bit of a floating-point value -* -* Inputs: -* val - the number -* -* Outputs: -* 0 if positive, non-zero if negative -* -**************************************************************** -* -__signbit start - - csubroutine (10:val),0 - - lda val+8 - and #$8000 - sta val - - creturn 2:val - end - -**************************************************************** -* -* int __fpcompare(long double x, long double y, short mask); -* -* Compare two floating-point values, not signaling invalid -* if they are unordered. -* -* Inputs: -* x,y - values to compare -* mask - mask of bits as returned in X register from FCMP -* -* Outputs: -* 1 if x and y have one of the relations specified by mask -* 0 otherwise -* -**************************************************************** -* -__fpcompare start - - csubroutine (10:x,10:y,2:mask),0 - - tdc - clc - adc #x - pea 0 - pha - tdc - clc - adc #y - pea 0 - pha - FCMPX - txa - and mask - beq lb1 - lda #1 -lb1 sta mask - - creturn 2:mask - end - -**************************************************************** -* -* double acosh(double x); -* -* Returns the inverse hyperbolic cosine of x. -* -**************************************************************** -* -acosh start -acoshf entry -acoshl entry - using MathCommon2 - - csubroutine (10:x),0 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - lda x y = sqrt(x-1) - sta y - lda x+2 - sta y+2 - lda x+4 - sta y+4 - lda x+6 - sta y+6 - lda x+8 - sta y+8 - ph4 #one - ph4 #y - FSUBI - ph4 #y - FSQRTX - - lda x t1 = sqrt(x+1) - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - lda x+8 - sta t1+8 - ph4 #one - ph4 #t1 - FADDI - ph4 #t1 - FSQRTX - - ph4 #y t1 = ln(1+y*(y+t1)) - ph4 #t1 - FADDX - ph4 #y - ph4 #t1 - FMULX - ph4 #t1 - FLN1X - - lda t1+8 if t1 = +inf - cmp #32767 - bne ret - lda t1+6 - asl a - ora t1+4 - ora t1+2 - ora t1 - bne ret - - pea 0 clear exceptions - FSETENV - lda x t1 = ln(x) + ln(2) - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - lda x+8 - sta t1+8 - ph4 #t1 - FLNX - ph4 #ln2 - ph4 #t1 - FADDX - -ret FPROCEXIT restore env & raise any new exceptions - plb - creturn 10:t1 return t1 - -y ds 10 temporary variable -one dc i'1' constants -ln2 dc e'0.69314718055994530942' - end - -**************************************************************** -* -* double asinh(double x); -* -* Returns the inverse hyperbolic sine of x. -* -**************************************************************** -* -asinh start -asinhf entry -asinhl entry - using MathCommon2 - - csubroutine (10:x),0 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - pei x+8 save sign of x - asl x+8 x = abs(x) - lsr x+8 - - lda x t1 = y = z = x - sta y - sta z - sta t1 - lda x+2 - sta y+2 - sta z+2 - sta t1+2 - lda x+4 - sta y+4 - sta z+4 - sta t1+4 - lda x+6 - sta y+6 - sta z+6 - sta t1+6 - lda x+8 - sta y+8 - sta z+8 - sta t1+8 - - lda x if value is zero (or typical inf) - ora x+2 - ora x+4 - ora x+6 - beq skipcalc return the input value - - lda x+8 else if x is very small - cmp #-33+16383 - bge calc - pea INEXACT raise "inexact" exception - FSETXCP -skipcalc brl setsign return the input value - -calc cmp #16383/2+16383 else if x is very large (or nan) - blt notbig - ph4 #z z = ln(x) + ln(2) - FLNX - ph4 #ln2 - ph4 #z - FADDX - brl setsign else - -notbig pea -2 t1 = 1 / (t1 * t1) - ph4 #t1 - FXPWRI - - ph4 #one t1 = 1 + t1 - ph4 #t1 - FADDI - - ph4 #t1 t1 = sqrt(t1) - FSQRTX - - pea -1 y = 1 / y - ph4 #y - FXPWRI - - ph4 #y t1 = t1 + y - ph4 #t1 - FADDX - - ph4 #t1 z = z / t1 - ph4 #z - FDIVX - - tdc z = z + x - clc - adc #x - pea 0 - pha - ph4 #z - FADDX - - ph4 #z z = ln(1+z) - FLN1X - -setsign asl z+8 sign of z = original sign of x - pla - asl a - ror z+8 - - FPROCEXIT restore env & raise any new exceptions - plb - creturn 10:z return z - -y ds 10 temporary variables -z ds 10 -one dc i'1' constants -ln2 dc e'0.69314718055994530942' - end - -**************************************************************** -* -* double atanh(double x); -* -* Returns the inverse hyperbolic tangent of x. -* -**************************************************************** -* -atanh start -atanhf entry -atanhl entry - using MathCommon2 - - csubroutine (10:x),0 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - pei x+8 save sign of x - asl x+8 x = abs(x) - lsr x+8 - - lda x t1 = x - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - lda x+8 - sta t1+8 - - lda x+8 if x is very small - cmp #-33+16383 - bge calc - lda x if value is not zero - ora x+2 - ora x+4 - ora x+6 - beq skipcalc - pea INEXACT raise "inexact" exception - FSETXCP -skipcalc bra setsign skip next steps (return input value) - -calc ph4 #one x = x - 1 - tdc - clc - adc #x - pea 0 - pha - FSUBI - - tdc t1 = t1 / x - clc - adc #x - pea 0 - pha - ph4 #t1 - FDIVX - - lda t1+8 if t1 is inf/nan - asl a - cmp #32767*2 - beq setsign skip next steps (so atanh(1) = +inf) - - ph4 #minustwo t1 = t1 * -2 - ph4 #t1 - FMULI - - ph4 #t1 t1 = ln(1+t1) - FLN1X - - ph4 #minustwo t1 = t1 / -2 - ph4 #t1 - FDIVI - -setsign asl t1+8 sign of t1 = original sign of x - pla - asl a - ror t1+8 - - FPROCEXIT restore env & raise any new exceptions - plb - creturn 10:t1 return t1 - -one dc i'1' constants -minustwo dc i'-2' - end - -**************************************************************** -* -* double cbrt(double x); -* -* Returns x^(1/3) (the cube root of x). -* -**************************************************************** -* -cbrt start -cbrtf entry -cbrtl entry - using MathCommon2 -scale equ 1 - - csubroutine (10:x),2 - - phb - phk - plb - - stz scale scale by 0 by default (for inf/nan) - - lda x+8 - pha save original sign - and #$7FFF - sta x+8 force sign to + - cmp #32767 skip scaling for inf/nan - beq do_calc - - ldx x+6 if number is denormalized - bmi div_exp - bne normaliz - ldx x+4 - bne normaliz - ldx x+2 - bne normaliz - ldx x - beq div_exp - -normaliz dec a normalize it and adjust exponent - asl x - rol x+2 - rol x+4 - rol x+6 - bpl normaliz - -div_exp pha calculate exponent/3 - pha - pha - pea 3 - _SDivide - pla a = quotient - plx x = remainder - cpx #2 adjust remainder of 2 to -1 - bne setscale - ldx #-1 - inc a - -setscale sec calculate amount to scale result by - sbc #16383/3 - sta scale - txa use remainder as exponent for calc. - clc - adc #16383 -do_calc sta t1+8 - - lda x place mantissa in work area - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - - ph4 #onethird compute val^(1/3) - ph4 #t1 - FXPWRY - - clc apply scaling - lda t1+8 - adc scale - sta t1+8 - - asl t1+8 set sign of result to orig. sign of x - pla - asl a - ror t1+8 - - plb - creturn 10:t1 return t1 - -onethird dc e'0.33333333333333333333' - end - -**************************************************************** -* -* double copysign(double x, double y); -* -* Returns a value with the magnitude of x and the sign of y. -* -**************************************************************** -* -copysign start -copysignf entry -copysignl entry - using MathCommon2 - - phb place x in a work area... - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - asl a ...with the sign bit shifted off - sta t1+8 - - pla remove y - pla - pla - pla - pla - asl a get sign bit of y - ror t1+8 give return value that sign - - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double erf(double x); -* -* Returns the error function of x. -* -* double erfc(double x); -* -* Returns the complementary error function of x, 1 - erf(x). -* -* This implementation is based on W. J. Cody's article -* "Rational Chebyshev Approximations for the Error Function." -* -**************************************************************** -* -erf start -erff entry -erfl entry - using MathCommon2 -x_offset equ 9 stack offset of x (in most of the code) - - clc - bra lb1 - -erfc entry -erfcf entry -erfcl entry - - sec - -lb1 phb save & set data bank - phk - plb - - pha make space for saved SANE environment - - lda x_offset-2+8,s - ror a save erf/erfc flag (high bit) - pha and original sign (next bit) - - rol a t1 := |x| - rol a - lsr a - sta t1+8 - lda x_offset+6,s - sta t1+6 - lda x_offset+4,s - sta t1+4 - lda x_offset+2,s - sta t1+2 - lda x_offset,s - sta t1 - - tsc save env & set to default - clc - adc #3 - pea 0 - pha - FPROCENTRY -; -; Computation using approximation of erf, for small enough x values -; - ph4 #threshold if |x| <= 0.5 then - ph4 #t1 - FCMPS - jmi use_erfc - - ph4 #t1 t1 := x^2 - ph4 #t1 - FMULX - - ph4 #y y := P1(t1) - pea 4 - ph4 #P1 - jsl poly - - ph4 #z z := Q1(t1) - pea 4 - ph4 #Q1 - jsl poly - - ph4 #z y := y/z - ph4 #y - FDIVX - - tsc y := x * y - clc - adc #x_offset - pea 0 - pha - ph4 #y - FMULX - - pla - jpl clearxcp if computing erfc then - - ph4 #one y := y - 1 - ph4 #y - FSUBX - brl flipsign y := -y -; -; Computation using approximations of erfc, for larger x values -; -use_erfc ph4 #four else - ph4 #t1 - FCMPI - jmi big_erfc if |x| <= 4 then - - ph4 #y y := P2(t1) - pea 8 - ph4 #P2 - jsl poly - - ph4 #z z := Q2(t1) - pea 8 - ph4 #Q2 - jsl poly - - ph4 #z y := y/z - ph4 #y - FDIVX - - ph4 #t1 t1 := e^(-x^2) - ph4 #t1 - FMULX - lda t1+8 - eor #$8000 - sta t1+8 - ph4 #t1 - FEXPX - - ph4 #t1 y := t1 * y - ph4 #y - FMULX - - brl end_erfc else (if |x| > 4 or NAN) - -big_erfc pea -2 t1 := 1 / x^2 - ph4 #t1 - FXPWRI - - ph4 #y y := P3(t1) - pea 5 - ph4 #P3 - jsl poly - - ph4 #z z := Q3(t1) - pea 5 - ph4 #Q3 - jsl poly - - ph4 #z y := y/z - ph4 #y - FDIVX - - ph4 #t1 y := t1 * y - ph4 #y - FMULX - - ph4 #one_over_sqrt_pi y := 1/sqrt(pi) + y - ph4 #y - FADDX - - lda x_offset+8,s y := y / |x| - and #$7fff - sta x_offset+8,s - tsc - clc - adc #x_offset - ldx #0 - phx (push operands of below calls) - pha - phx - pha - phx - pha - phx - pha - phx - pha - ph4 #y - FDIVX - - FMULX y := e^(-x^2) * y - lda x_offset+8+8,s - eor #$8000 - sta x_offset+8+8,s - FEXPX - ph4 #y - FMULX - -end_erfc pla - bpl erf_from_erfc if computing erfc then - - ldx #$1300 (set allowed exception mask) - asl a - bpl rstr_env if x < 0 - - ph4 #two y := y - 2 - ph4 #y - FSUBI - bra flipsign y := -y - -erf_from_erfc anop - pha - ph4 #one if computing erf then - ph4 #y - FSUBX y := y - 1 - - pla - asl a - bmi clearxcp if x > 0 then - -flipsign lda y+8 y := -y - eor #$8000 - sta y+8 - -clearxcp ldx #$1100 ignore overflow, div-by-zero -rstr_env stx z (& underflow unless doing erfc for x>.5) - FGETENV - txa - and z - ora 1,s - sta 1,s - FSETENV unless computing erfc for x > 4 - - pla clean up stack - sta 9,s - pla - sta 9,s - tsc - clc - adc #6 - tcs - plb - ldx #^y return a pointer to the result - lda #y - rtl - -threshold dc f'0.5' threshold for computing erf or erfc - -; constants -two dc i2'2' -four dc i2'4' -one_over_sqrt_pi dc e'0.564189583547756286924' -; coefficients for erf calculation, |x| <= .5 -P1 dc e'1.857777061846031526730e-1' - dc e'3.161123743870565596947e+0' - dc e'1.138641541510501556495e+2' - dc e'3.774852376853020208137e+2' - dc e'3.209377589138469472562e+3' -one anop -Q1 dc e'1.0' - dc e'2.360129095234412093499e+1' - dc e'2.440246379344441733056e+2' - dc e'1.282616526077372275645e+3' - dc e'2.844236833439170622273e+3' -; coefficients for erfc calculation, .46875 <= x <= 4 -P2 dc e'2.15311535474403846343e-8' - dc e'5.64188496988670089180e-1' - dc e'8.88314979438837594118e+0' - dc e'6.61191906371416294775e+1' - dc e'2.98635138197400131132e+2' - dc e'8.81952221241769090411e+2' - dc e'1.71204761263407058314e+3' - dc e'2.05107837782607146532e+3' - dc e'1.23033935479799725272e+3' -Q2 dc e'1.0' - dc e'1.57449261107098347253e+1' - dc e'1.17693950891312499305e+2' - dc e'5.37181101862009857509e+2' - dc e'1.62138957456669018874e+3' - dc e'3.29079923573345962678e+3' - dc e'4.36261909014324715820e+3' - dc e'3.43936767414372163696e+3' - dc e'1.23033935480374942043e+3' -; coefficients for erfc calculation, x >= 4 -P3 dc e'-1.63153871373020978498e-2' - dc e'-3.05326634961232344035e-1' - dc e'-3.60344899949804439429e-1' - dc e'-1.25781726111229246204e-1' - dc e'-1.60837851487422766278e-2' - dc e'-6.58749161529837803157e-4' -Q3 dc e'1.0' - dc e'2.56852019228982242072e+0' - dc e'1.87295284992346047209e+0' - dc e'5.27905102951428412248e-1' - dc e'6.05183413124413191178e-2' - dc e'2.33520497626869185443e-3' -; temporaries / return values -y ds 10 -z ds 10 - end - -**************************************************************** -* -* double exp2(double x); -* -* Returns 2^x. -* -**************************************************************** -* -exp2 start -exp2f entry -exp2l entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FEXP2X - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double expm1(double x); -* -* Returns e^x - 1. -* -**************************************************************** -* -expm1 start -expm1f entry -expm1l entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FEXP1X - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double fdim(double x, double y); -* -* Returns x - y if x > y, or +0 if x <= y. -* -**************************************************************** -* -fdim start -fdimf entry -fdiml entry - using MathCommon2 - - phb - phk - plb - - tsc compare x and y - clc - adc #5 - pea 0 - pha - adc #10 - pea 0 - pha - FCMPX - bmi x_le_y - beq x_le_y - - tsc if x > y (or unordered) - clc - adc #5+10 - pea 0 - pha - sbc #10-1 (carry is clear) - pea 0 - pha - FSUBX x = x - y - lda 5,s t1 = x - sta t1 - lda 5+2,s - sta t1+2 - lda 5+4,s - sta t1+4 - lda 5+6,s - sta t1+6 - lda 5+8,s - sta t1+8 - bra ret else - -x_le_y stz t1 t1 = +0.0 - stz t1+2 - stz t1+4 - stz t1+6 - stz t1+8 - -ret plx clean up stack - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double fma(double x, double y, double z); -* -* Compute (x * y) + z, rounded only once at the end. -* -**************************************************************** -* -fma start -fmaf entry -fmal entry - using MathCommon2 -mant1 equ 1 mantissa of value 1 -exp1 equ mant1+16 exponent of value 1 -sign1 equ exp1+4 sign of value 1 (high bit) -mant2 equ sign1+2 mantissa of value 2 -exp2 equ mant2+16 exponent of value 2 -sign2 equ exp2+4 sign of value 2 (low bit) -expdiff equ sign2+2 difference between exponents -extra equ expdiff+4 extra bits (guard, round, sticky) -xcps equ extra+2 floating-point exceptions - - csubroutine (10:x,10:y,10:z),54 - - stz extra - lda x if x or y is NAN, INF, or 0 then - ora x+2 - ora x+4 - ora x+6 - beq nanInf0 return (x * y) + z computed with SANE - lda x+8 - asl a - cmp #32767*2 - beq nanInf0 - - lda y - ora y+2 - ora y+4 - ora y+6 - beq nanInf0 - lda y+8 - asl a - cmp #32767*2 - beq nanInf0 - - lda z+8 else if z is INF or NAN then - asl a - cmp #32767*2 - beq x_plus_z return x + z computed with SANE - inc extra - lda z else if z is 0 then - ora z+2 - ora z+4 - ora z+6 return x * y computed with SANE - bne compute else compute fma(x,y,z) ourselves - -; -; Compute with SANE if any operands are NAN/INF/0 -; -nanInf0 tdc if in first or third case above then - clc - adc #y - pea 0 - pha - adc #x-y - pea 0 - pha - FMULX x = x * y - -x_plus_z ldy extra if in first or second case above then - bne return_x - - tdc - clc - adc #z - phy - pha - adc #x-z - phy - pha - FADDX x = x + z - -return_x lda x copy result to t1 - sta >t1 - lda x+2 - sta >t1+2 - lda x+4 - sta >t1+4 - lda x+6 - sta >t1+6 - lda x+8 - sta >t1+8 - brl ret return result - -; -; Compute it ourselves if all operands are finite and non-zero -; -compute stz xcps no exceptions so far - lda x copy mantissa of x to mant1 - sta mant1 - lda x+2 - sta mant1+2 - lda x+4 - sta mant1+4 - lda x+6 - sta mant1+6 - stz mant1+8 - stz mant1+10 - stz mant1+12 - stz mant1+14 - - ldy #64 multiply mantissas (64 x 64 to 128-bit) -ml1 lda mant1 - lsr a - bcc ml2 - clc add multiplicand to partial product - lda mant1+8 - adc y - sta mant1+8 - lda mant1+10 - adc y+2 - sta mant1+10 - lda mant1+12 - adc y+4 - sta mant1+12 - lda mant1+14 - adc y+6 - sta mant1+14 -ml2 ror mant1+14 shift the interim result - ror mant1+12 - ror mant1+10 - ror mant1+8 - ror mant1+6 - ror mant1+4 - ror mant1+2 - ror mant1 - dey loop until done - bne ml1 - - lda x+8 calculate exponent - asl a - sta exp1 - lda y+8 - asl a - clc - adc exp1 - ror a - sta exp1 - stz exp1+2 - add4 exp1,#-16383+1 - - lda mant1+14 normalize calculated value - bmi getsign1 -norm1_lp dec4 exp1 - asl mant1 - rol mant1+2 - rol mant1+4 - rol mant1+6 - rol mant1+8 - rol mant1+10 - rol mant1+12 - rol mant1+14 - bpl norm1_lp - -getsign1 lda x+8 get sign of x*y - eor y+8 - sta sign1 - - lda z+8 get sign of z - sta sign2 - - and #$7fff copy exponent of z to exp2 - sta exp2 - stz exp2+2 - - stz mant2 copy mantissa of z to mant2 - stz mant2+2 - stz mant2+4 - stz mant2+6 - lda z - sta mant2+8 - lda z+2 - sta mant2+10 - lda z+4 - sta mant2+12 - lda z+6 - sta mant2+14 - - bmi exp_cmp normalize z value -norm2_lp dec4 exp2 - asl mant2+8 (low mantissa bits stay 0) - rol mant2+10 - rol mant2+12 - rol mant2+14 - bpl norm2_lp - -exp_cmp cmp4 exp1,exp2 if exp1 < exp2 - bge do_align - jsr exchange exchange value 1 and value 2 - -; at this point, exp1 >= exp2 -do_align stz extra initially extra bits are 0 - sub4 exp1,exp2,expdiff expdiff = exp1 - exp2 - cmpl expdiff,#65+1 if expdiff > 65 then - blt aligntst - stz mant2 zero out mant2 - stz mant2+2 - stz mant2+4 - stz mant2+6 - stz mant2+8 - stz mant2+10 - stz mant2+12 - stz mant2+14 - inc extra but set the sticky bit for rounding - bra addorsub else - -align_lp dec4 expdiff - lsr mant2+14 shift mant2 until it is aligned - ror mant2+12 - ror mant2+10 - ror mant2+8 - ror mant2+6 - ror mant2+4 - ror mant2+2 - ror mant2 - ror extra - bcc aligntst maintain sticky bit - lda #$0001 - tsb extra -aligntst lda expdiff - ora expdiff+2 - bne align_lp - -addorsub lda sign1 if signs of x*y and z are the same then - eor sign2 - bmi subtract - - clc mant1 = mant1 + mant2 - ldx #-16 -addLoop lda mant1+16,x - adc mant2+16,x - sta mant1+16,x - inx - inx - bmi addLoop - bcc add_done if there is carry out - ror mant1+14 rotate carry back into result - ror mant1+12 - ror mant1+10 - ror mant1+8 - ror mant1+6 - ror mant1+4 - ror mant1+2 - ror mant1 - ror extra - bcc inc_exp maintain sticky bit - lda #$0001 - tsb extra -inc_exp inc4 exp1 increment exponent -add_done bra xtrabits else - -subtract ldx #14 if mant1 < mant2 then -subCmpLp lda mant1,x (note: only occurs if mant2 was - cmp mant2,x not shifted, so extra is 0) - bne sub_cmp - dex - dex - bpl subCmpLp -sub_cmp bge do_sub - jsr exchange exchange mant2 and mant1 - -do_sub sec mant1 = mant1 - mant2 (including extra) - lda #0 - sbc extra - sta extra - ldx #-16 -subLoop lda mant1+16,x - sbc mant2+16,x - sta mant1+16,x - inx - inx - bmi subLoop - ora mant1 if result (including extra bits) is 0 then - ora mant1+2 - ora mant1+4 - ora mant1+6 - ora mant1+8 - ora mant1+10 - ora mant1+12 - ora extra - bne subalign - stz exp1 set exponent to 0 - stz sign1 set sign to + - FGETENV if rounding direction is downward then - txa - bpl savezero - asl a - bmi savezero - dec sign1 set sign to - -savezero brl do_save skip to return -subalign lda mant1+14 - bmi xtrabits normalize after subtraction, if needed -subAl_lp dec4 exp1 - asl extra - rol mant1 - rol mant1+2 - rol mant1+4 - rol mant1+6 - rol mant1+8 - rol mant1+10 - rol mant1+12 - rol mant1+14 -subAlNeg bpl subAl_lp - -xtrabits lda mant1 consolidate extra bits (into mant1+6) - ora mant1+2 - ora mant1+4 - ora extra - beq denorm - lda #$0001 - tsb mant1+6 - -denorm lda #INEXACT assume INEXACT is just INEXACT - bra denormCk while exponent is too small -denormLp inc4 exp1 increment exponent - lsr mant1+14 shift mantissa right - ror mant1+12 - ror mant1+10 - ror mant1+8 - ror mant1+6 - bcc denorm2 maintain sticky bit - lda #$0001 - tsb mant1+6 -denorm2 lda #UNDERFLOW+INEXACT flag that INEXACT also implies UNDERFLOW -denormCk ldy exp1+2 - bmi denormLp - - ldy mant1+6 if there are extra bits then - beq saveval - tsb xcps set inexact (+ maybe underflow) exception - FGETENV get rounding direction - txa - asl a - bcs roundDn0 - bmi roundUp if rounding to nearest then - lda mant1+6 if first extra bit is 0 - bpl saveval do not round - asl a else if remaining extra bits are non-zero - bne do_round - lda mant1+8 or low-order bit of result is 1 then - lsr a - bcc saveval - bra do_round apply rounding - -roundUp lda sign1 if rounding upward then - bmi saveval if positive then - bra do_round apply rounding - -roundDn0 bmi saveval if rounding downward then - lda sign1 if negative then - bpl saveval apply rounding - -do_round inc mant1+8 (perform the rounding, if needed) - bne saveval - inc mant1+10 - bne saveval - inc mant1+12 - bne saveval - inc mant1+14 - bne saveval - sec - ror mant1+14 - ror mant1+12 - ror mant1+10 - ror mant1+8 - inc4 exp1 - -saveval lda exp1+2 if value is too large to represent then - bne save_inf - lda exp1 - cmp #32766+1 - blt do_save -save_inf lda #32767 set it to infinity - sta exp1 - stz mant1+8 - stz mant1+10 - stz mant1+12 - stz mant1+14 - lda #OVERFLOW+INEXACT set overflow and inexact exceptions - tsb xcps -do_save lda mant1+8 generate result - sta >t1 - lda mant1+10 - sta >t1+2 - lda mant1+12 - sta >t1+4 - lda mant1+14 - sta >t1+6 - lda exp1 - asl a - asl sign1 - ror a - sta >t1+8 - - lda xcps if there were exceptions then - beq ret - pha set them in SANE environment - FSETXCP - -ret creturn 10:t1 return t1 - -; local subroutine - exchange value 1 and value 2 -; Note: requires mant1/exp1/sign1 and mant2/exp2/sign2 to be in order -exchange ldx #16+4+2-2 -xchgLp lda mant1,x - ldy mant2,x - sta mant2,x - sty mant1,x - dex - dex - bpl xchgLp - rts - end - -**************************************************************** -* -* double fmax(double x, double y); -* -* Returns the maximum numeric value of x or y. -* If one is a NaN, returns the other. -* -**************************************************************** -* -fmax start -fmaxf entry -fmaxl entry - using MathCommon2 - - phb - phk - plb - phd - - tsc set up direct page - clc - adc #7 - tcd - - pea 0 compare x and y - pha - clc - adc #10 - pea 0 - pha - FCMPX - - bmi use_y if x < y, return y - bvs use_x if x >= y, return x - beq use_x - - pea 0 if x,y are unordered - phd - FCLASSX - txa - and #$00FE - cmp #$00FC if x is not a nan, return x - beq use_y else return y - -use_x ldx #0 - bra copyit - -use_y ldx #10 - -copyit lda 0,x copy result to t1 - sta t1 - lda 2,x - sta t1+2 - lda 4,x - sta t1+4 - lda 6,x - sta t1+6 - lda 8,x - sta t1+8 - - pld clean up stack - plx - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double fmin(double x, double y); -* -* Returns the minimum numeric value of x or y. -* If one is a NaN, returns the other. -* -**************************************************************** -* -fmin start -fminf entry -fminl entry - using MathCommon2 - - phb - phk - plb - phd - - tsc set up direct page - clc - adc #7 - tcd - - pea 0 compare x and y - pha - clc - adc #10 - pea 0 - pha - FCMPX - - bmi use_x if x < y, return x - bvs use_y if x >= y, return y - beq use_y - - pea 0 if x,y are unordered - phd - FCLASSX - txa - and #$00FE - cmp #$00FC if x is not a nan, return x - beq use_y else return y - -use_x ldx #0 - bra copyit - -use_y ldx #10 - -copyit lda 0,x copy result to t1 - sta t1 - lda 2,x - sta t1+2 - lda 4,x - sta t1+4 - lda 6,x - sta t1+6 - lda 8,x - sta t1+8 - - pld clean up stack - plx - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double hypot(double x, double y); -* -* Returns the square root of x^2 + y^2, without undue overflow -* or underflow. -* -**************************************************************** -* -hypot start -hypotf entry -hypotl entry - using MathCommon2 -scale equ 1 scaling factor - - csubroutine (10:x,10:y),2 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - stz scale no scaling by default - - asl x+8 x = abs(x) - lsr x+8 - asl y+8 y = abs(y) - lsr y+8 - - tdc if x < y - clc - adc #x - pea 0 - pha - adc #y-x - pea 0 - pha - FCMPX - bpl sorted - - ldx #8 exchange x and y -xchgloop lda x,x - ldy y,x - sta y,x - sty x,x - dex - dex - bpl xchgloop -sorted anop at this point, 0 <= y <= x (if ordered) - - lda x+8 if x or y is nan or inf - ldy y+8 - cpy #32767 - beq naninf - cmp #32767 - beq naninf skip exponent manipulation - - cmp #8190+16383+1 if exponent of x > 8190 - blt chksmall - sec scale x and y down by 2^8300 - sbc #8300 - sta x+8 - lda #8300 - sta scale - lda y+8 - sec - sbc #8300 - sta y+8 - bpl compute - stz y (zero out y if needed) - stz y+2 - stz y+4 - stz y+6 - stz y+8 - bra compute - -chksmall cmp #-8100+16383 else if exponent of x < -8100 - bge compute - clc scale x and y up by 2^8300 - adc #8300 - sta x+8 - lda y+8 - clc - adc #8300 - sta y+8 - lda #-8300 - sta scale - -compute tdc x = x*x - clc - adc #x - pea 0 - pha - pea 0 - pha - FMULX - - tdc y = y*y - clc - adc #y - pea 0 - pha - pea 0 - pha - FMULX - -naninf anop (we skip to here if x or y is nan/inf) - lda x copy x to t1 - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - lda x+8 - sta t1+8 - - tdc t1 = x*x + y*y - clc - adc #y - pea 0 - pha - ph4 #t1 - FADDX - - ph4 #t1 t1 = sqrt(t1) - FSQRTX - - lda scale if scaling is needed - beq done - pha do it - ph4 #t1 - FSCALBX - -done FPROCEXIT restore env - plb - creturn 10:t1 return t1 - end - -**************************************************************** -* -* int ilogb(double x); -* -* Returns the binary exponent of x (a signed integer value), -* treating denormalized numbers as if they were normalized. -* Handles inf/nan/0 cases specially. -* -**************************************************************** -* -ilogb start -ilogbf entry -ilogbl entry - - csubroutine (10:x),0 - - tdc check for special cases - clc - adc #x - pea 0 - pha - FCLASSX - ldy #$7FFF - txa - and #$FF - cmp #$FE if x is INF - beq special return INT_MAX - lsr a - beq do_logb if x is 0 or NAN - iny return INT_MIN -special sty x - bra ret - -do_logb tdc compute logb(x) - clc - adc #x - pea 0 - pha - FLOGBX - - tdc convert to integer - clc - adc #x - pea 0 - pha - pea 0 - pha - FX2I - -ret creturn 2:x return it - rtl - end - -**************************************************************** -* -* double lgamma(double x); -* -* Computes the natural logarithm of the absolute value of the -* gamma function of x. -* -* The rational approximations used near 1 and 2 are from -* W. J. Cody and K. E. Hillstrom, "Chebyshev Approximations -* for the Natural Logarithm of the Gamma Function". -* For other positive x values, a computation based on -* Stirling's formula is used. -* -**************************************************************** -* -lgamma start -lgammaf entry -lgammal entry - using MathCommon2 - - csubroutine (10:x),0 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - stz reflect assume no reflection - -; For x < 0, lgamma(x) = ln(pi) - ln(abs(x*sin(pi*x))) - lgamma(-x) - lda x+8 if x is negative then - jpl positive - - lda x if x != -0 and x != -inf then - ora x+2 - ora x+4 - ora x+6 - jeq abs_x - - inc reflect flag to do reflection - - ldx #8 neg_adj := x -lp1 lda x,X - sta neg_adj,X - dex - dex - bpl lp1 - - ph4 #two neg_adj := neg_adj REM 2 - ph4 #neg_adj - FREMI - - ph4 #pi neg_adj := sin(neg_adj*pi) - ph4 #neg_adj - FMULX - ph4 #neg_adj - FSINX - - tdc neg_adj := neg_adj * x - clc - adc #x - pea 0 - pha - ph4 #neg_adj - FMULX - - asl neg_adj+8 neg_adj := abs(neg_adj) - lsr neg_adj+8 - - ph4 #neg_adj neg_adj := ln(neg_adj) - FLNX - - ph4 #ln_pi neg_adj := neg_adj - ln(pi) - ph4 #neg_adj - FSUBX - - lda neg_adj+8 if neg_adj is inf/nan then - asl a - cmp #32767*2 - bne abs_x - - ldx #8 z := neg_adj -lp2 lda neg_adj,X - sta z,X - dex - dex - bpl lp2 - - brl ret_negz return -z - -abs_x asl x+8 x := abs(x) - lsr x+8 - -positive ldx #8 t1 := x -lp3 lda x,X - sta t1,X - dex - dex - bpl lp3 - - ph4 #rat2_hi if x is near 1 or 2 then - ph4 #t1 - FCMPS - jvc not_near_1_or_2 - -; Rational approximation for x near 2 - ph4 #rat2_lo if x is near 2 then - ph4 #t1 - FCMPS - bvs small - - ph4 #z z := P2(t1) - pea 7 - ph4 #P2 - jsl poly - - ph4 #y y := Q2(t1) - pea 7 - ph4 #Q2 - jsl poly - - ph4 #two N := 2 - bra finish_rational_approx - -; Rational approximation for x near 1 -small ph4 #rat1_hi else if x is near 1 then - ph4 #t1 - FCMPS - jvc not_near_1_or_2 - - ph4 #rat1_lo - ph4 #t1 - FCMPS - bvs not_near_1_or_2 - - ph4 #z z := P1(t1) - pea 7 - ph4 #P1 - jsl poly - - ph4 #y y := Q1(t1) - pea 7 - ph4 #Q1 - jsl poly - - ph4 #one N := 1 - -finish_rational_approx anop - ph4 #t1 t1 := x - N - FSUBI - - ph4 #y z := z/y - ph4 #z - FDIVX - - ph4 #t1 z := z * t1 - ph4 #z - FMULX - brl chk_reflect goto chk_reflect - -; Calculate based on Stirling's formula for other x values -not_near_1_or_2 anop - stz scaled assume no scaling used - - stz z z := 1.0 - stz z+2 - stz z+4 - lda #$8000 - sta z+6 - lda #16383 - sta z+8 - -chk_big ph4 #cutoff while x < 9.41796875 do - tdc - clc - adc #x - pea 0 - pha - FCMPS - bvc stirling - - inc scaled flag that scaling was used - - tdc z := z * x - clc - adc #x - pea 0 - pha - ph4 #z - FMULX - - ph4 #one x := x + 1 - tdc - clc - adc #x - pea 0 - pha - FADDI - bra chk_big - -stirling ldx #8 t1 := x -lp4 lda x,X y := x - sta t1,X - sta y,X - dex - dex - bpl lp4 - - ph4 #t1 t1 := ln(t1) - FLNX - - ph4 #one_half y := y - 1/2 - ph4 #y - FSUBS - - ph4 #t1 y := t1 * y - ph4 #y - FMULX - - lda y+8 if y is not inf/nan then - asl a - cmp #32767*2 - beq stirling_poly - - tdc y := y - x - clc - adc #x - pea 0 - pha - ph4 #y - FSUBX - - ph4 #half_ln_2pi y := y + 0.5*ln(2pi) - ph4 #y - FADDX - - lda scaled if scaled then - beq stirling_poly - - ph4 #z z := ln(z) - FLNX - - ph4 #z y := y - z - ph4 #y - FSUBX - -stirling_poly anop - ldx #8 t1 := x -lp5 lda x,X - sta t1,X - dex - dex - bpl lp5 - - pea -2 t1 := 1/t1^2 - ph4 #t1 - FXPWRI - - ph4 #z z := P_stirling(t1) - pea 8 - ph4 #P_stirling - jsl poly - - tdc z := z / x - clc - adc #x - pea 0 - pha - ph4 #z - FDIVX - - ph4 #y z := y + z - ph4 #z - FADDX - -chk_reflect anop - lda reflect if reflection is needed then - beq done - - ph4 #neg_adj z := z + neg_adj - ph4 #z - FADDX - -ret_negz lda z+8 z := -z - eor #$8000 - sta z+8 - -done FPROCEXIT restore env & raise any new exceptions - plb - creturn 10:z return a pointer to the result - -; constants -cutoff dc f'9.41796875' above this, just use Stirling's formula -;low/high limits for use of rational approximations near 1 and 2 -rat1_lo dc f'0.97265625' -rat1_hi dc f'1.02880859375' -rat2_lo dc f'1.76220703125' -rat2_hi dc f'2.208984375' - -one_half dc f'0.5' -half_ln_2pi dc e'0.9189385332046727417803297' -one dc i2'1' -two dc i2'2' -ln_pi dc e'1.14472988584940017414342735' -pi dc e'3.1415926535897932384626433' - -; coefficients for rational approximation, 0.5 <= x <= 1.5 -P1 dc e'4.120843185847770031e00' - dc e'8.568982062831317339e01' - dc e'2.431752435244210223e02' - dc e'-2.617218583856145190e02' - dc e'-9.222613728801521582e02' - dc e'-5.176383498023217924e02' - dc e'-7.741064071332953034e01' - dc e'-2.208843997216182306e00' -Q1 dc e'1.000000000000000000e00' - dc e'4.564677187585907957e01' - dc e'3.778372484823942081e02' - dc e'9.513235976797059772e02' - dc e'8.460755362020782006e02' - dc e'2.623083470269460180e02' - dc e'2.443519662506311704e01' - dc e'4.097792921092615065e-01' -; coefficients for rational approximation, 1.5 <= x <= 4.0 -P2 dc e'5.1550576176408171704e00' - dc e'3.7751067979721702241e02' - dc e'5.2689832559149812458e03' - dc e'1.9553605540630449846e04' - dc e'1.2043173809871640151e04' - dc e'-2.0648294205325283281e04' - dc e'-1.5086302287667250272e04' - dc e'-1.5138318341150667785e03' -Q2 dc e'1.0000000000000000000e00' - dc e'1.2890931890129576873e02' - dc e'3.0399030414394398824e03' - dc e'2.2029562144156636889e04' - dc e'5.7120255396025029854e04' - dc e'5.2622863838411992470e04' - dc e'1.4402090371700852304e04' - dc e'6.9832741405735102159e02' -; coefficients for Stirling series approximation -P_stirling anop - dc e'0.1796443723688305731649385' - dc e'-0.02955065359477124183006536' - dc e'0.006410256410256410256410256' - dc e'-0.001917526917526917526917527' - dc e'0.0008417508417508417508417508' - dc e'-0.0005952380952380952380952381' - dc e'0.0007936507936507936507936508' - dc e'-0.002777777777777777777777778' - dc e'0.08333333333333333333333333' - -; temporaries/return values -y ds 10 -z ds 10 - -neg_adj ds 10 adjustment for negative x -reflect ds 2 reflection flag -scaled ds 2 scaling flag - end - -**************************************************************** -* -* long long llrint(double x); -* -* Rounds x to an integer using current rounding direction -* and returns it as a long long (if representable). -* -* Note: This avoids calling FX2C on negative numbers, -* because it is buggy for certain values. -* -**************************************************************** -* -llrint start -llrintf entry -llrintl entry -retptr equ 1 - - csubroutine (10:x),4 - stx retptr - stz retptr+2 - - tdc - clc - adc #x - pea 0 push src address for fcpxx - pha - pea llmin|-16 push dst address for fcpxx - pea llmin - pea 0 push operand address for frintx - pha - FRINTX round - FCPXX compare with LLONG_MIN - bne convert - - lda #$8000 if it is LLONG_MIN, use that value - ldy #6 - sta [retptr],y - asl a - dey - dey - sta [retptr],y - dey - dey - sta [retptr],y - sta [retptr] - bra done otherwise - -convert pei x+8 save sign of x - asl x+8 x = abs(x) - lsr x+8 - tdc - clc - adc #x - pea 0 push src address for fx2c - pha - pei retptr+2 push dst address for fx2c - pei retptr - FX2C convert x - - pla if x was negative - bpl done - sec - lda #0 negate result - sbc [retptr] - sta [retptr] - ldy #2 - lda #0 - sbc [retptr],y - sta [retptr],y - iny - iny - lda #0 - sbc [retptr],y - sta [retptr],y - iny - iny - lda #0 - sbc [retptr],y - sta [retptr],y - -done creturn - -llmin dc e'-9223372036854775808' - end - -**************************************************************** -* -* long long llround(double x); -* -* Rounds x to the nearest integer, rounding halfway cases away -* from 0, and returns it as a long long (if representable). -* -* Note: This avoids calling FX2C on negative numbers, -* because it is buggy for certain values. -* -**************************************************************** -* -llround start -llroundf entry -llroundl entry -retptr equ 1 - - csubroutine (10:x),4 - stx retptr - stz retptr+2 - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - tdc if x == LLONG_MIN - clc - adc #x - pea 0 - pha - ph4 #llmin - FCMPX - beq retllmin return LLONG_MIN - tdc else if x == LLONG_MIN+0.5 - clc - adc #x - pea 0 - pha - ph4 #llminp05 - FCPXX - bne convert - - pea INEXACT raise "inexact" exception - FSETXCP -retllmin lda #$8000 return LLONG_MIN - ldy #6 - sta [retptr],y - asl a - dey - dey - sta [retptr],y - dey - dey - sta [retptr],y - sta [retptr] - brl ret else - -convert pei x+8 save sign of x - asl x+8 x = abs(x) - lsr x+8 - - tdc round to integer - clc - adc #x - pea 0 - pha - pei retptr+2 - pei retptr - FX2C - - pea INEXACT - FTESTXCP if there was no inexact exception - beq chk_neg we're done: x was an integer/nan/inf - - FGETENV else - txa - ora #TOWARDZERO*$4000 round toward zero - pha - FSETENV - - ph4 #onehalf x = x + 0.5 (rounded toward 0) - tdc - clc - adc #x - pea 0 - pha - FADDS - tdc round to integer - clc - adc #x - pea 0 - pha - pei retptr+2 - pei retptr - FX2C - -chk_neg pla if x was negative - bpl ret - sec - lda #0 negate result - sbc [retptr] - sta [retptr] - ldy #2 - lda #0 - sbc [retptr],y - sta [retptr],y - iny - iny - lda #0 - sbc [retptr],y - sta [retptr],y - iny - iny - lda #0 - sbc [retptr],y - sta [retptr],y - -ret FPROCEXIT restore env & raise any new exceptions - creturn - -llmin dc e'-9223372036854775808' -llminp05 dc e'-9223372036854775807.5' -onehalf dc f'0.5' - end - -**************************************************************** -* -* double log1p(double x); -* -* Returns ln(1+x). -* -**************************************************************** -* -log1p start -log1pf entry -log1pl entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FLN1X - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double log2(double x); -* -* Returns log2(x) (the base-2 logarithm of x). -* -**************************************************************** -* -log2 start -log2f entry -log2l entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FLOG2X - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double logb(double x); -* -* Returns the binary exponent of x (a signed integer value), -* treating denormalized numbers as if they were normalized. -* -**************************************************************** -* -logb start -logbf entry -logbl entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FLOGBX - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* long lrint(double x); -* -* Rounds x to an integer using current rounding direction -* and returns it as a long (if representable). -* -* Note: This avoids calling FX2L or FX2C on negative numbers, -* because they are buggy for certain values. -* -**************************************************************** -* -lrint start -lrintf entry -lrintl entry - - csubroutine (10:x),0 - - pei x+8 save sign of x - - tdc - clc - adc #x - pea 0 - pha - pea 0 - pha - pea 0 - pha - FRINTX round x to integer - asl x+8 x = abs(x) - lsr x+8 - FX2C convert to comp - - lda x+4 if x is out of range of long - ora x+6 - bne flag_inv - cmpl x,#$80000000 - blt chk_neg - bne flag_inv - lda 1,s - bmi chk_neg -flag_inv pea INVALID raise "invalid" exception - FSETXCP - -chk_neg pla if x was negative - bpl ret - sub4 #0,x,x negate result - -ret creturn 4:x return it - rtl - end - -**************************************************************** -* -* long lround(double x); -* -* Rounds x to the nearest integer, rounding halfway cases -* away from 0, and returns it as a long (if representable). -* -* Note: This avoids calling FX2L or FX2C on negative numbers, -* because they are buggy for certain values. -* -**************************************************************** -* -lround start -lroundf entry -lroundl entry -result equ 1 result value - - csubroutine (10:x),8 - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - - pei x+8 save sign of x - asl x+8 x = abs(x) - lsr x+8 - - tdc round to integer with default rounding - clc - adc #x - pea 0 - pha - adc #result-x - pea 0 - pha - FX2C - - pea INEXACT - FTESTXCP if there was no inexact exception - beq chkrange we are done: x was an integer/nan/inf - - FGETENV - txa - ora #TOWARDZERO*$4000 set rounding direction to "toward zero" - pha - FSETENV - - ph4 #onehalf x = x + 0.5 (rounded toward 0) - tdc - clc - adc #x - pea 0 - pha - FADDS - tdc round to integer - clc - adc #x - pea 0 - pha - adc #result-x - pea 0 - pha - FX2C - -chkrange lda result+4 if x is out of range of long - ora result+6 - bne flag_inv - cmpl result,#$80000000 - blt chk_neg - bne flag_inv - lda 1,s - bmi chk_neg -flag_inv pea INVALID raise "invalid" exception - FSETXCP - -chk_neg pla if x was negative - bpl ret - sub4 #0,result,result negate result - -ret FPROCEXIT restore env & raise any new exceptions - creturn 4:result return the result - -onehalf dc f'0.5' - end - -**************************************************************** -* -* float modff(float x, float *iptr); -* -* Splits x into integer and fractional parts. Returns the -* fractional part and stores integer part as a float in *iptr. -* -**************************************************************** -* -modff start - using MathCommon2 - - csubroutine (10:x,4:iptr),0 - - phb - phk - plb - - lda x copy x to t1 - sta t1 - lda x+2 - sta t1+2 - lda x+4 - sta t1+4 - lda x+6 - sta t1+6 - lda x+8 - sta t1+8 - - asl a check for infinity or nan - cmp #32767|1 - bne finite - lda x+6 - asl a - ora x+4 - ora x+2 - ora x - bne storeint if value is nan, return it as-is - stz t1 if value is +-inf, fractional part is 0 - stz t1+2 - stz t1+4 - stz t1+6 - stz t1+8 - bra storeint - -finite tdc truncate x to an integer - clc - adc #x - pea 0 - pha - FTINTX - - tdc t1 := t1 - x - clc - adc #x - pea 0 - pha - ph4 #t1 - FSUBX - -storeint tdc copy x to *iptr, converting to float - clc - adc #x - pea 0 - pha - pei iptr+2 - pei iptr - FX2S - -copysign asl t1+8 copy sign of x to t1 - asl x+8 - ror t1+8 - - plb - creturn 10:t1 return t1 (fractional part) - end - -**************************************************************** -* -* long double modfl(long double x, long double *iptr); -* -* Splits x into integer and fractional parts. Returns the -* fractional part and stores the integer part in *iptr. -* -**************************************************************** -* -modfl start - using MathCommon2 - - csubroutine (10:x,4:iptr),0 - - phb - phk - plb - - lda x copy x to *iptr and t1 - sta [iptr] - sta t1 - ldy #2 - lda x+2 - sta [iptr],y - sta t1+2 - iny - iny - lda x+4 - sta [iptr],y - sta t1+4 - iny - iny - lda x+6 - sta [iptr],y - sta t1+6 - iny - iny - lda x+8 - sta [iptr],y - sta t1+8 - - asl a check for infinity or nan - cmp #32767|1 - bne finite - lda x+6 - asl a - ora x+4 - ora x+2 - ora x - bne ret if value is nan, return it as-is - stz t1 if value is +-inf, fractional part is 0 - stz t1+2 - stz t1+4 - stz t1+6 - stz t1+8 - bra copysign - -finite pei iptr+2 if value is finite - pei iptr - FTINTX truncate *iptr to an integer - - pei iptr+2 t1 := t1 - *iptr - pei iptr - ph4 #t1 - FSUBX - -copysign asl t1+8 copy sign of x to t1 - asl x+8 - ror t1+8 - -ret plb - creturn 10:t1 return t1 (fractional part) - end - -**************************************************************** -* -* double nan(const char *tagp); -* -* Returns a quiet NaN, with NaN code determined by the -* argument string. -* -**************************************************************** -* -nan start -nanf entry -nanl entry - using MathCommon2 - - csubroutine (4:tagp) - - phb - phk - plb - - stz t1+6 initial code is 0 - -loop lda [tagp] do - and #$00FF get next character - beq loopdone if end of string, break - cmp #'0' - blt no_code - cmp #'9'+1 - bge no_code if not a digit, treat as no code - and #$000F - asl t1+6 code = code*10 + digit - clc - adc t1+6 - asl t1+6 - asl t1+6 - clc - adc t1+6 - sta t1+6 - inc4 tagp tagp++ - bra loop while true - -no_code stz t1+6 if no code specified, default to 0 - -loopdone lda t1+6 - and #$00FF use low 8 bits as NaN code - bne codeok if code is 0 - lda #21 use NANZERO -codeok ora #$4000 set high bit of f for quiet NaN - sta t1+6 - - lda #32767 e=32767 for NaN - sta t1+8 - stz t1+4 set rest of fraction field to 0 - stz t1+2 - stz t1 - - plb - creturn 10:t1 return a pointer to the result - end - -**************************************************************** -* -* double nearbyint(double x); -* -* Rounds x to an integer using current rounding direction, -* never raising the "inexact" exception. -* -**************************************************************** -* -nearbyint start -nearbyintf entry -nearbyintl entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - FGETENV save environment - phx - ph4 #t1 compute the value - FRINTX - FSETENV restore environment - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double nextafter(double x, double y); -* -* Returns next representable value (in double format) -* after x in the direction of y. Returns y if x equals y. -* -**************************************************************** -* -nextafter start - using MathCommon2 - - tsc x = (double) x - clc - adc #4 - pea 0 - pha - pea 0 - pha - FX2D - lda 4,s save low bits of x - sta 4+8,s - - tsc y = (double) y - clc - adc #4+10 - pea 0 - pha - pea 0 - pha - FX2D - - tsc push address of y - clc - adc #4+10 - pea 0 - pha - sbc #10-1 push address of x - pea 0 - pha - FNEXTD x = nextafter x toward y - - tsc store x (as extended) in t1 - clc - adc #4 - pea 0 - pha - ph4 #t1 - FD2X - - phb - lda 4+8+1,s if original x might be 0 then - bne ret - tsc - clc - adc #4+10+1 - pea 0 - pha - ph4 #t1 - FCPXD - bne ret if t1 == y then - phk - plb - asl t1+8 sign of t1 = sign of y - lda 4+10+1+6,s - asl a - ror t1+8 - -ret plx move return address - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* float nextafterf(float x, float y); -* -* Returns next representable value (in float format) -* after x in the direction of y. Returns y if x equals y. -* -**************************************************************** -* -nextafterf start - using MathCommon2 - - tsc x = (float) x - clc - adc #4 - pea 0 - pha - pea 0 - pha - FX2S - lda 4,s save low bits of x - sta 4+8,s - - tsc y = (float) y - clc - adc #4+10 - pea 0 - pha - pea 0 - pha - FX2S - - tsc push address of y - clc - adc #4+10 - pea 0 - pha - sbc #10-1 push address of x - pea 0 - pha - FNEXTS x = nextafter x toward y - - tsc store x (as extended) in t1 - clc - adc #4 - pea 0 - pha - ph4 #t1 - FS2X - - phb - lda 4+8+1,s if original x might be 0 then - bne ret - tsc - clc - adc #4+10+1 - pea 0 - pha - ph4 #t1 - FCPXS - bne ret if t1 == y then - phk - plb - asl t1+8 sign of t1 = sign of y - lda 4+10+1+2,s - asl a - ror t1+8 - -ret plx move return address - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* long double nextafterl(long double x, long double y); -* long double nexttowardl(long double x, long double y); -* -* Returns next representable value (in extended format) -* after x in the direction of y. Returns y if x equals y. -* -**************************************************************** -* -nextafterl start -nexttowardl entry - using MathCommon2 - - tsc push address of x - clc - adc #4 - pea 0 - pha - adc #10 push address of y - pea 0 - pha - FCPXX - bne getnext if x == y then - tsc - clc - adc #4+10 return y - bra storeval else - -getnext tsc push address of y - clc - adc #4+10 - pea 0 - pha - sbc #10-1 push address of x - pea 0 - pha - FNEXTX x = nextafter x toward y - - tsc return x - clc - adc #4 -storeval pea 0 store return value to t1 - pha - ph4 #t1 - FX2X - - phb move return address - plx - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double nexttoward(double x, long double y); -* -* Returns next representable value (in double format) -* after x in the direction of y. Returns y if x equals y. -* -**************************************************************** -* -nexttoward start - using MathCommon2 - - tsc x = (double) x - clc - adc #4 - pea 0 - pha - pea 0 - pha - FX2D - - tsc push address of x - clc - adc #4 - pea 0 - pha - adc #10 push address of y - pea 0 - pha - FCPXD compare x and y - - bvs x_gt_y - bmi x_lt_y - beq x_eq_y - - tsc x,y unordered case: do nextafter(x,y) - clc - adc #4+10 - pea 0 - pha - pea 0 - pha - pea 0 - pha - FX2D - bra getnext - -x_gt_y ph4 #minusinf x > y case: do nextafter(x,-inf) - bra getnext - -x_lt_y ph4 #plusinf x < y case: do nextafter(x,+inf) - bra getnext - -x_eq_y phb - phk - plb - lda 4+10+1,s x == y case: return y - sta t1 - lda 4+10+1+2,s - sta t1+2 - lda 4+10+1+4,s - sta t1+4 - lda 4+10+1+6,s - sta t1+6 - lda 4+10+1+8,s - sta t1+8 - bra ret - -getnext tsc compute nextafter(x,...) - clc - adc #4+4 - pea 0 - pha - FNEXTD - - tsc store x (as extended) in t1 - clc - adc #4 - pea 0 - pha - ph4 #t1 - FD2X - - phb move return address -ret plx - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - -plusinf dc d'+inf' -minusinf dc d'-inf' - end - -**************************************************************** -* -* float nexttowardf(float x, long double y); -* -* Returns next representable value (in float format) -* after x in the direction of y. Returns y if x equals y. -* -**************************************************************** -* -nexttowardf start - using MathCommon2 - - tsc x = (double) x - clc - adc #4 - pea 0 - pha - pea 0 - pha - FX2S - - tsc push address of x - clc - adc #4 - pea 0 - pha - adc #10 push address of y - pea 0 - pha - FCPXS compare x and y - - bvs x_gt_y - bmi x_lt_y - beq x_eq_y - - tsc x,y unordered case: do nextafter(x,y) - clc - adc #4+10 - pea 0 - pha - pea 0 - pha - pea 0 - pha - FX2S - bra getnext - -x_gt_y ph4 #minusinf x > y case: do nextafter(x,-inf) - bra getnext - -x_lt_y ph4 #plusinf x < y case: do nextafter(x,+inf) - bra getnext - -x_eq_y phb - phk - plb - lda 4+10+1,s x == y case: return y - sta t1 - lda 4+10+1+2,s - sta t1+2 - lda 4+10+1+4,s - sta t1+4 - lda 4+10+1+6,s - sta t1+6 - lda 4+10+1+8,s - sta t1+8 - bra ret - -getnext tsc compute nextafter(x,...) - clc - adc #4+4 - pea 0 - pha - FNEXTS - - tsc store x (as extended) in t1 - clc - adc #4 - pea 0 - pha - ph4 #t1 - FS2X - - phb move return address -ret plx - ply - tsc - clc - adc #20 - tcs - phy - phx - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - -plusinf dc f'+inf' -minusinf dc f'-inf' - end - -**************************************************************** -* -* poly: evaluate a polynomial -* -* Evaluates sum from i=0 to n of K_i * x^i -* -* Inputs: -* coeffs: array of coefficients, K_n down to K_0 -* n: degree of polynomial -* result: pointer to location for result -* t1: x value -* -* Note: The coeffs array is assumed not to cross banks. -* -**************************************************************** -* -poly private - using MathCommon2 - - csubroutine (4:coeffs,2:n,4:result),0 - - ldy #8 val := K_n -loop1 lda [coeffs],y - sta [result],y - dey - dey - bpl loop1 - -loop2 lda coeffs for i := n-1 downto 0 - clc - adc #10 - sta coeffs - - ph4 #t1 val := val * x - ph4 INT_MAX - blt notbig - pea 32767 scale by INT_MAX - pea 0 - bra adjust_n -notbig cmp4 n,#-32768 else if n < INT_MIN - bge notsmall - pea -32768+64 scale by INT_MIN - pea -1 - -adjust_n sec if n is out of range of int - lda n subtract scale factor from n - sbc 3,s - sta n - lda n+2 - sbc 1,s - sta n+2 - pla - bra do_scalb else -notsmall pei n scale by n - stz n remaining amount to scale by is 0 - stz n+2 - -do_scalb ph4 #t1 scale the number - FSCALBX - - lda n if no more scaling to do - ora n+2 - beq done we are done - - ph4 #t1 else if value is nan/inf/zero - FCLASSX - txa - and #$FE - bne done stop: more scaling would not change it - brl loop else scale by remaining amount - -done plb - creturn 10:t1 return a pointer to the result - end - -**************************************************************** -* -* double scalbn(double x, int n); -* -* Returns x * 2^n. -* -**************************************************************** -* -scalbn start -scalbnf entry -scalbnl entry - using MathCommon2 - - phb place x in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - - pla get n - phy - phx - - pha compute the value - ph4 #t1 - FSCALBX - - plb - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* double tgamma(double x); -* -* Computes the gamma function of x. -* -**************************************************************** -* -tgamma start -tgammaf entry -tgammal entry - using MathCommon2 - - csubroutine (10:x),0 - - phb - phk - plb - - pha save env & set to default - tsc - inc a - pea 0 - pha - FPROCENTRY - -* For x < 0.5, use Gamma(x) = pi / (sin(x * pi) * Gamma(1 - x)) - stz reflect - ph4 #one_half if x < 0.5 then - tdc - clc - adc #x - pea 0 - pha - FCMPS - bvc lb1 - inc reflect - - ldx #8 orig_x := x -lp0 lda x,x - sta fracpart,x - sta orig_x,x - dex - dex - bpl lp0 - - ph4 #two fracpart := x REM 2 - ph4 #fracpart - FREMI - - ph4 #one x := x-1 - tdc - clc - adc #x - pea 0 - pha - FSUBX - - lda x+8 - eor #$8000 - sta x+8 - -* For 0 <= x <= 9.375, use the identity Gamma(x) = Gamma(x+1)/x -lb1 ldy #8 denom := 1 -lp1 lda one,y - sta denom,y - dey - dey - bpl lp1 - - ph4 #cutoff - ph4 #z - FS2X - - tdc z := 10.375 - x - clc - adc #x - pea 0 - pha - ph4 #z - FSUBX - - lda z+8 if z < 0 (or NAN) - jmi stirling just use Stirling series approx. - cmp #$7fff - jeq stirling - - ph4 #z truncate z to integer - FTINTX - - ph4 #z x := x + z - tdc - clc - adc #x - pea 0 - pha - FADDX - - tdc y := x - clc - adc #x - pea 0 - pha - ph4 #y - FX2X - - ph4 #z repeat z times : - ph4 #z - FX2I - -lp2 dec z - bmi stirling - - ph4 #one y := y - 1 - ph4 #y - FSUBX - - ph4 #y denom := denom * y - ph4 #denom - FMULX - - bra lp2 - -* For x >= 9.375, calculate Gamma(x) using a Stirling series approximation -stirling lda x t1 := x - sta y y := x - sta z z := x - sta t1 - lda x+2 - sta y+2 - sta z+2 - sta t1+2 - lda x+4 - sta y+4 - sta z+4 - sta t1+4 - lda x+6 - sta y+6 - sta z+6 - sta t1+6 - lda x+8 - sta y+8 - sta z+8 - sta t1+8 - - ph4 #e z := x/e - ph4 #z - FDIVX - ph4 #one_half y := x - 1/2 - ph4 #y - FSUBS - ph4 #y z := (x/e)^(x-1/2) - ph4 #z - FXPWRY - - pea -1 t1 := 1/x - ph4 #t1 - FXPWRI - ph4 #y y := P(t1) - pea 17 - ph4 #P - jsl poly - - ph4 #y z := z * y * sqrt(2*pi/e) - ph4 #z - FMULX - ph4 #sqrt_2pi_over_e - ph4 #z - FMULX - -* Adjust result as necessary for small or negative x values - ph4 #denom z := z / denom - ph4 #z (for cases where initial x was small) - FDIVX - - lda reflect if doing reflection - jeq done - - ph4 #pi fracpart := sin(x*pi) - ph4 #fracpart - FMULX - ph4 #fracpart - FSINX - - ph4 #fracpart if sin(x*pi)=0 (i.e. x was an integer) - FCLASSX - txa - inc a - and #$00ff - bne lb2 - - asl fracpart+8 take sign from original x (for +-0) - asl orig_x+8 - ror fracpart+8 - - lda orig_x if original x was not 0 - ora orig_x+2 - ora orig_x+4 - ora orig_x+6 - beq lb2 - - lda #32767 force NAN result - sta z+8 - sta z+6 - - pea $0100 raise "invalid" exception (only) - FSETENV - -lb2 ph4 #z z := pi / (fracpart * z) - ph4 #fracpart - FMULX - - ph4 #pi - ph4 #z - FX2X - - ph4 #fracpart - ph4 #z - FDIVX - -done FPROCEXIT restore env & raise any new exceptions - plb - creturn 10:z return a pointer to the result - -cutoff dc f'10.375' cutoff for Stirling approximation (+1) - -one_half dc f'0.5' -two dc i2'2' -e dc e'2.7182818284590452353602874713526624977572' -pi dc e'3.1415926535897932384626433' -sqrt_2pi_over_e dc e'1.520346901066280805611940146754975627' - -P anop Stirling series constants - dc e'+1.79540117061234856108e-01' - dc e'-2.48174360026499773092e-03' - dc e'-2.95278809456991205054e-02' - dc e'+5.40164767892604515180e-04' - dc e'+6.40336283380806979482e-03' - dc e'-1.62516262783915816899e-04' - dc e'-1.91443849856547752650e-03' - dc e'+7.20489541602001055909e-05' - dc e'+8.39498720672087279993e-04' - dc e'-5.17179090826059219337e-05' - dc e'-5.92166437353693882865e-04' - dc e'+6.97281375836585777429e-05' - dc e'+7.84039221720066627474e-04' - dc e'-2.29472093621399176955e-04' - dc e'-2.68132716049382716049e-03' - dc e'+3.47222222222222222222e-03' - dc e'+8.33333333333333333333e-02' -one dc e'+1.00000000000000000000e+00' - -y ds 10 -z ds 10 -denom ds 10 -fracpart ds 10 - -reflect ds 2 flag: do reflection? -orig_x ds 10 original x value - end - -**************************************************************** -* -* double trunc(double x); -* -* Truncates x to an integer (discarding fractional part). -* -**************************************************************** -* -trunc start -truncf entry -truncl entry - using MathCommon2 - - phb place the number in a work area - plx - ply - phk - plb - pla - sta t1 - pla - sta t1+2 - pla - sta t1+4 - pla - sta t1+6 - pla - sta t1+8 - phy - phx - plb - - ph4 #t1 compute the value - FTINTX - - ldx #^t1 return a pointer to the result - lda #t1 - rtl - end - -**************************************************************** -* -* float and long double versions of functions in SysFloat -* -**************************************************************** -* -acosf start -acosl entry - jml acos - end - -asinf start -asinl entry - jml asin - end - -atanf start -atanl entry - jml atan - end - -atan2f start -atan2l entry - jml atan2 - end - -ceilf start -ceill entry - jml ceil - end - -cosf start -cosl entry - jml cos - end - -coshf start -coshl entry - jml cosh - end - -expf start -expl entry - jml exp - end - -fabsf start -fabsl entry - jml fabs - end - -floorf start -floorl entry - jml floor - end - -fmodf start -fmodl entry - jml fmod - end - -frexpf start -frexpl entry - jml frexp - end - -ldexpf start -ldexpl entry - jml ldexp - end - -logf start -logl entry - jml log - end - -log10f start -log10l entry - jml log10 - end - -powf start -powl entry - jml pow - end - -sinf start -sinl entry - jml sin - end - -sinhf start -sinhl entry - jml sinh - end - -sqrtf start -sqrtl entry - jml sqrt - end - -tanf start -tanl entry - jml tan - end - -tanhf start -tanhl entry - jml tanh - end diff --git a/math2.macros b/math2.macros deleted file mode 100644 index 9ab9e35..0000000 --- a/math2.macros +++ /dev/null @@ -1,755 +0,0 @@ - macro -&l ph4 &n1 -&l anop - aif "&n1"="*",.f - lclc &c -&c amid &n1,1,1 - aif "&c"="#",.d - aif s:longa=1,.a - rep #%00100000 -.a - aif "&c"<>"{",.b -&c amid &n1,l:&n1,1 - aif "&c"<>"}",.g -&n1 amid &n1,2,l:&n1-2 - ldy #2 - lda (&n1),y - pha - lda (&n1) - pha - ago .e -.b - aif "&c"<>"[",.c - ldy #2 - lda &n1,y - pha - lda &n1 - pha - ago .e -.c - aif "&c"<>"<",.c1 -&n1 amid &n1,2,l:&n1-1 - pei &n1+2 - pei &n1 - ago .e -.c1 - lda &n1+2 - pha - lda &n1 - pha - ago .e -.d -&n1 amid &n1,2,l:&n1-1 - pea +(&n1)|-16 - pea &n1 - ago .f -.e - aif s:longa=1,.f - sep #%00100000 -.f - mexit -.g - mnote "Missing closing '}'",16 - mend - MACRO -&lab csubroutine &parms,&work -&lab anop - aif c:&work,.a - lclc &work -&work setc 0 -.a - gbla &totallen - gbla &worklen -&worklen seta &work -&totallen seta 0 - aif c:&parms=0,.e - lclc &len - lclc &p - lcla &i -&i seta 1 -.b -&p setc &parms(&i) -&len amid &p,2,1 - aif "&len"=":",.c -&len amid &p,1,2 -&p amid &p,4,l:&p-3 - ago .d -.c -&len amid &p,1,1 -&p amid &p,3,l:&p-2 -.d -&p equ &totallen+4+&work -&totallen seta &totallen+&len -&i seta &i+1 - aif &i<=c:&parms,^b -.e - tsc - aif &work=0,.f - sec - sbc #&work - tcs -.f - phd - tcd - mend - MACRO -&lab creturn &r -&lab anop - lclc &len - aif c:&r,.a - lclc &r -&r setc 0 -&len setc 0 - ago .h -.a -&len amid &r,2,1 - aif "&len"=":",.b -&len amid &r,1,2 -&r amid &r,4,l:&r-3 - ago .c -.b -&len amid &r,1,1 -&r amid &r,3,l:&r-2 -.c - aif &len<>2,.d - ldy &r - ago .h -.d - aif &len<>4,.e - ldx &r+2 - ldy &r - ago .h -.e - aif &len<>10,.g - ldy #&r - ldx #^&r - ago .h -.g - mnote 'Not a valid return length',16 - mexit -.h - aif &totallen=0,.i - lda &worklen+2 - sta &worklen+&totallen+2 - lda &worklen+1 - sta &worklen+&totallen+1 -.i - pld - tsc - clc - adc #&worklen+&totallen - tcs - aif &len=0,.j - tya -.j - rtl - mend - macro -&l cmp4 &n1,&n2 - lclb &yistwo -&l ~setm - ~lda.h &n1 - ~op.h eor,&n2 - bpl ~a&SYSCNT - ~lda.h &n2 - ~op.h cmp,&n1 - bra ~b&SYSCNT -~a&SYSCNT ~lda.h &n1 - ~op.h cmp,&n2 - bne ~b&SYSCNT - ~lda &n1 - ~op cmp,&n2 -~b&SYSCNT anop - ~restm - mend - macro -&l ~lda &op - lclc &c -&c amid "&op",1,1 - aif "&c"<>"{",.b -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b -&l lda &op - mend - macro -&l ~lda.h &op -&l anop - lclc &c -&c amid "&op",1,1 - aif "&c"="[",.b - aif "&c"<>"{",.d -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b - aif &yistwo,.c -&yistwo setb 1 - ldy #2 -.c -&op setc "&op,y" - lda &op - mexit -.d - aif "&c"<>"#",.e -&op amid "&op",2,l:&op-1 -&op setc "#^&op" - lda &op - mexit -.e - lda 2+&op - mend - macro -&l ~op &opc,&op - lclc &c -&c amid "&op",1,1 - aif "&c"<>"{",.b -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b -&l &opc &op - mend - macro -&l ~op.h &opc,&op -&l anop - lclc &c -&c amid "&op",1,1 - aif "&c"="[",.b - aif "&c"<>"{",.d -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b - aif &yistwo,.c -&yistwo setb 1 - ldy #2 -.c -&op setc "&op,y" - &opc &op - mexit -.d - aif "&c"<>"#",.e -&op amid "&op",2,l:&op-1 -&op setc "#^&op" - &opc &op - mexit -.e - &opc 2+&op - mend - macro -&l ~restm -&l anop - aif (&~la+&~li)=2,.i - sep #32*(.not.&~la)+16*(.not.&~li) - aif &~la,.h - longa off -.h - aif &~li,.i - longi off -.i - mend - macro -&l ~setm -&l anop - aif c:&~la,.b - gblb &~la - gblb &~li -.b -&~la setb s:longa -&~li setb s:longi - aif s:longa.and.s:longi,.a - rep #32*(.not.&~la)+16*(.not.&~li) - longa on - longi on -.a - mend - macro -&l inc4 &a -&l ~setm - inc &a - bne ~&SYSCNT - inc 2+&a -~&SYSCNT ~restm - mend - macro -&l sub4 &m1,&m2,&m3 - lclb &yistwo - lclc &c -&l ~setm - aif c:&m3,.a -&c amid "&m2",1,1 - aif "&c"<>"#",.a -&c amid "&m1",1,1 - aif "&c"="{",.a - aif "&c"="[",.a -&c amid "&m2",2,l:&m2-1 - aif &c>=65536,.a - sec - ~lda &m1 - ~op sbc,&m2 - ~sta &m1 - bcs ~&SYSCNT - ~op.h dec,&m1 -~&SYSCNT anop - ago .c -.a - aif c:&m3,.b - lclc &m3 -&m3 setc &m1 -.b - sec - ~lda &m1 - ~op sbc,&m2 - ~sta &m3 - ~lda.h &m1 - ~op.h sbc,&m2 - ~sta.h &m3 -.c - ~restm - mend - macro -&l ~sta &op - lclc &c -&c amid "&op",1,1 - aif "&c"<>"{",.b -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b -&l sta &op - mend - macro -&l ~sta.h &op -&l anop - lclc &c -&c amid "&op",1,1 - aif "&c"="[",.b - aif "&c"<>"{",.d -&c amid "&op",l:&op,1 - aif "&c"="}",.a - mnote "Missing closing '}'",2 -&op setc &op} -.a -&op amid "&op",2,l:&op-2 -&op setc (&op) -.b - aif &yistwo,.c -&yistwo setb 1 - ldy #2 -.c -&op setc "&op,y" - sta &op - mexit -.d - sta 2+&op - mend - macro -&l cmpl &n1,&n2 - lclb &yistwo -&l ~setm - ~lda.h &n1 - ~op.h cmp,&n2 - bne ~a&SYSCNT - ~lda &n1 - ~op cmp,&n2 -~a&SYSCNT anop - ~restm - mend - macro -&l jmi &bp -&l bpl *+5 - brl &bp - mend - macro -&l jpl &bp -&l bmi *+5 - brl &bp - mend - macro -&l jeq &bp -&l bne *+5 - brl &bp - mend - MACRO -&LAB FCLASSS -&LAB PEA $021C - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCLASSD -&LAB PEA $011C - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCLASSX -&LAB PEA $001C - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2S -&LAB PEA $0210 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2D -&LAB PEA $0110 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCMPX -&LAB PEA $0008 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FEXP2X -&LAB PEA $000A - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FEXP1X -&LAB PEA $000C - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FLN1X -&LAB PEA $0004 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FLOG2X -&LAB PEA $0002 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FLOGBX -&LAB PEA $001A - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2I -&LAB PEA $0410 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FTINTX -&LAB PEA $0016 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FRINTX -&LAB PEA $0014 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FXPWRY -&LAB PEA $0012 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FREMX -&LAB PEA $000C - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSCALBX -&LAB PEA $0018 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSUBX -&LAB PEA $0002 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FGETENV -&LAB PEA $03 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSETENV -&LAB PEA $01 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2C -&LAB PEA $0510 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCPXX -&LAB PEA $0A - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FNEXTX -&LAB PEA $001E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FX2X -&LAB PEA $0010 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCPXD -&LAB PEA $010A - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FNEXTD -&LAB PEA $011E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FD2X -&LAB PEA $010E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FS2X -&LAB PEA $020E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FNEXTS -&LAB PEA $021E - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FCPXS -&LAB PEA $020A - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FPROCENTRY -&LAB PEA $0017 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FPROCEXIT -&LAB PEA $0019 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FTESTXCP -&LAB PEA $001B - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FADDS -&LAB PEA $0200 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSETXCP -&LAB PEA $0015 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FADDX -&LAB PEA $0000 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FADDI -&LAB PEA $0400 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSUBI -&LAB PEA $0402 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FMULX -&LAB PEA $0004 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSQRTX -&LAB PEA $0012 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FLNX -&LAB PEA $0000 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&lab _SDivide -&lab ldx #$0A0B - jsl $E10000 - MEND - MACRO -&LAB FMULI -&LAB PEA $0404 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FDIVI -&LAB PEA $0406 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FDIVX -&LAB PEA $0006 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FXPWRI -&LAB PEA $0010 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FCMPS -&LAB PEA $0208 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FEXPX -&LAB PEA $0008 - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FCMPI -&LAB PEA $0408 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSUBS -&LAB PEA $0202 - LDX #$090A - JSL $E10000 - MEND - MACRO -&LAB FSINX -&LAB PEA $001A - LDX #$0B0A - JSL $E10000 - MEND - MACRO -&LAB FREMI -&LAB PEA $040C - LDX #$090A - JSL $E10000 - MEND - macro -&l dec4 &a -&l ~setm - lda &a - bne ~&SYSCNT - dec 2+&a -~&SYSCNT dec &a - ~restm - mend - macro -&l add4 &m1,&m2,&m3 - lclb &yistwo - lclc &c -&l ~setm - aif c:&m3,.a -&c amid "&m2",1,1 - aif "&c"<>"#",.a -&c amid "&m1",1,1 - aif "&c"="{",.a - aif "&c"="[",.a -&c amid "&m2",2,l:&m2-1 - aif &c>=65536,.a - clc - ~lda &m1 - ~op adc,&m2 - ~sta &m1 - bcc ~&SYSCNT - ~op.h inc,&m1 -~&SYSCNT anop - ago .c -.a - aif c:&m3,.b - lclc &m3 -&m3 setc &m1 -.b - clc - ~lda &m1 - ~op adc,&m2 - ~sta &m3 - ~lda.h &m1 - ~op.h adc,&m2 - ~sta.h &m3 -.c - ~restm - mend - macro -&l jvc &bp -&l bvs *+5 - brl &bp - mend diff --git a/stdio.asm b/stdio.asm index 1ecfa4a..e499bb4 100644 --- a/stdio.asm +++ b/stdio.asm @@ -3828,289 +3828,6 @@ ug1 rtl string ds 4 end -**************************************************************** -* -* ~Format_a - format a floating-point number in hex format -* (lowercase output) -* ~Format_A - format a floating-point number in hex format -* (uppercase output) -* -* Inputs: -* ~altForm - always include decimal point? -* ~fieldWidth - output field width -* ~paddChar - padd character -* ~leftJustify - left justify the output? -* ~precision - precision of output -* ~precisionSpecified - was the precision specified? -* -**************************************************************** -* -~Format_a private - using ~printfCommon -argp equ 7 argument pointer -; -; Set the "or" value; this is used to set the case of character results -; - lda #$20 - sta ~orVal - bra in1 - -~Format_A entry - stz ~orVal -; -; Check for infinities or nans -; -in1 ldy #8 load exponent/sign word - lda [argp],Y - asl a - tax - eor #32767*2 - bne sn1 if number is an infinity or NaN - lda #' ' do not use '0' padding - sta ~paddChar - lda ~orVal if doing %A format - bne in2 - brl ~Format_E format like %E -in2 brl ~Format_e else format like %e -; -; Determine sign -; -sn1 stz ~sgn assume sign is positive - bcc ex1 if sign of number is negative - lda #'-' set sign character to '-' - sta ~sign - dec ~sgn flag that sign is negative -; -; Get exponent -; -ex1 txa get exponent field - lsr a (clears carry) - sbc #16383-1 compute unbiased exponent - sta ~exp save it -; -; Get significand -; -sg1 ldy #6 store the significand in ~sig -sg2 lda [argp],Y - sta ~sig,Y - dey - dey - bpl sg2 - ora ~sig+2 if significand is zero then - ora ~sig+4 - ora ~sig+6 - bne pc1 - lda #3 set exponent so it will print as 0 - sta ~exp -; -; Determine precision -; -pc1 lda ~precisionSpecified if the precision was not specified then - bne rd0 - lda #15 use a precision of 15 - sta ~precision -; -; Do rounding -; -rd0 lda ~precision if precision < 15 - cmp #15 - jge pd1 - - stz ~sig+8 make sure bit above significand is zero - inc a shift significand (precision+1)*4 bits left - asl a - asl a - pha - tay -rd1 ldx #0 - clc -rd1a rol ~sig,X - inx - inx - txa - eor #16 - bne rd1a - dey - bne rd1 - - lda ~sig consolidate extra bits - ora ~sig+2 - ora ~sig+4 - beq rd2 - lda #1 - tsb ~sig+6 - -rd2 lda ~sig+6 if there are extra non-zero bits then - beq rdW - FGETENV get rounding direction - txa - asl a - bcs roundDn0 - bmi roundUp if rounding to nearest then -roundNr lda ~sig+6 if first extra bit is 0 - bpl rdW do not round - asl a else if remaining extra bits are non-zero - bne do_round - lda ~sig+8 or low-order bit of result is 1 then - lsr a - bcc rdW - bra do_round apply rounding - -roundUp lda ~sgn if rounding upward then - bmi rdW if positive then - bra do_round apply rounding - -roundDn0 bmi rdW if rounding downward then -roundDn lda ~sgn if negative then - bpl rdW apply rounding - -do_round ldx #8 (perform the rounding, if needed) -rdV inc ~sig,X - bne rdW - inx - inx - cpx #14+1 - blt rdV - -rdW ply shift significand (precision+1)*4 bits right -rdX ldx #14 - clc -rdXa ror ~sig,X - dex - dex - bpl rdXa - dey - bne rdX - - lsr ~sig+8 handle carry out from rounding - bcc pd1 - ldx #6 -rdYa ror ~sig,X - dex - dex - bpl rdYa - inc ~exp -; -; Compute amount of padding -; -pd1 lda ~fieldWidth subtract off precision from field width - sec - sbc ~precision - sec subtract off minimal extra chars - sbc #6 - sta ~fieldWidth - lda ~sign if there is a sign character then - beq pd2 - dec ~fieldWidth decrement field width -pd2 lda ~precision if precision != 0 or # flag used then - ora ~altForm - beq pd2a - sta ~altForm flag this - dec ~fieldWidth decrement field width -pd2a lda ~exp get exponent - bpl pd3 compute absolute value of exponent - eor #$FFFF - inc a -pd3 cmp #10 if |exponent| >= 10 then - blt pd4 - dec ~fieldWidth decrement field width - cmp #100 if |exponent| >= 100 then - blt pd4 - dec ~fieldWidth decrement field width - cmp #1000 if |exponent| >= 1000 then - blt pd4 - dec ~fieldWidth decrement field width - cmp #10000 if |exponent| >= 10000 then - blt pd4 - dec ~fieldWidth decrement field width -pd4 lda ~paddChar if we are not padding with zeros then - cmp #'0' - beq pn1 - jsr ~RightJustify handle right justification -; -; Print the number -; -pn1 lda ~sign if there is a sign character then - beq pn2 - pha print it - jsl ~putchar -pn2 pea '0' print hex prefix - jsl ~putchar - lda #'X' - ora ~orVal - pha - jsl ~putchar - jsr ~ZeroPad pad with '0's if needed - -pn5 lda #0 print the digits - ldy #4 -pn6 asl ~sig - rol ~sig+2 - rol ~sig+4 - rol ~sig+6 - rol a - dey - bne pn6 -; clc (already clear) - adc #'0' - cmp #'9'+1 - blt pn7 - adc #6 - ora ~orVal -pn7 pha - jsl ~putchar - lda ~altForm print '.' after first digit if needed - beq pn8 - ph2 #'.' - jsl ~putchar - stz ~altForm -pn8 dec ~precision - bpl pn5 -; -; Print exponent -; - lda #'P' print 'P' or 'p' exponent prefix - ora ~orVal - pha - jsl ~putchar - - lda ~exp adjust exponent to reflect 4 bits - dec a in integer part (before '.') - dec a - dec a - pha push exponent - bmi pe1 print '+' if exponent is positive - ph2 #'+' - jsl ~putchar -pe1 ph4 #~str push the string addr - ph2 #6 push the string buffer length - ph2 #1 do a signed conversion - _Int2Dec convert exponent to string - ldx #0 print the exponent -pe2 lda ~str,x - and #$00FF - cmp #' '+1 - blt pe3 - phx - pha - jsl ~putchar - plx -pe3 inx - cpx #6 - blt pe2 -; -; Remove the number from the argument list -; - lda argp remove the parameter - adc #10-1 (carry is set) - sta argp -; -; Handle left justification -; - brl ~LeftJustify handle left justification - end - - **************************************************************** * * ~Format_c - format a '%' character diff --git a/stdio.macros b/stdio.macros index c95ac18..14ab942 100644 --- a/stdio.macros +++ b/stdio.macros @@ -974,14 +974,3 @@ phd tcd mend - MACRO -&LAB FGETENV -&LAB PEA $03 - LDX #$090A - JSL $E10000 - MEND - macro -&l jge &bp -&l blt *+5 - brl &bp - mend