diff --git a/Makefile b/Makefile index 6914cb9..50a7f52 100644 --- a/Makefile +++ b/Makefile @@ -10,8 +10,7 @@ LIB = orca SRCS = cc.asm ctype.asm orca.asm signal2.c stdlib.asm string.asm \ - time.asm toolglue.asm vars.asm int64.asm fenv.asm fpextra.asm \ - math2.asm locale.asm uchar.asm + time.asm toolglue.asm vars.asm int64.asm locale.asm uchar.asm buildall .PHONY: build assert.o 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 4bebb88..0000000 --- a/math2.asm +++ /dev/null @@ -1,3890 +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 - -**************************************************************** -* -* 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 7c3e1d9..0000000 --- a/math2.macros +++ /dev/null @@ -1,750 +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 diff --git a/orca.asm b/orca.asm index 41e28fc..265280c 100644 --- a/orca.asm +++ b/orca.asm @@ -53,7 +53,7 @@ lb1 rtl * enddesk start - jmp ~ENDDESK + brl ~ENDDESK end **************************************************************** @@ -64,7 +64,7 @@ enddesk start * endgraph start - jmp ~ENDGRAPH + brl ~ENDGRAPH end **************************************************************** @@ -119,7 +119,7 @@ id dc 8c' ',i1'0' * startdesk start - jmp ~STARTDESK + brl ~STARTDESK end **************************************************************** @@ -130,7 +130,7 @@ startdesk start * startgraph start - jmp ~STARTGRAPH + brl ~STARTGRAPH end **************************************************************** diff --git a/stdio.asm b/stdio.asm index 1b2917f..e499bb4 100644 --- a/stdio.asm +++ b/stdio.asm @@ -1322,11 +1322,10 @@ fprintf start sta stream+2 phy restore return address/data bank phx + ldx stream plb - lda >stream+2 verify that stream exists - pha - lda >stream - pha + pha verify that stream exists + phx jsl ~VerifyStream bcc lb1 lda #EIO @@ -3327,11 +3326,10 @@ vfprintf start sta stream+2 phy restore return address/data bank phx + ldx stream plb - lda >stream+2 verify that stream exists - pha - lda >stream - pha + pha verify that stream exists + phx jsl ~VerifyStream bcc lb1 lda #EIO @@ -3881,40 +3879,33 @@ argp equ 7 argument pointer ; ; For signed numbers, if the value is negative, use the sign flag ; - lda ~isLongLong handle long long values - beq sn0 - ldy #6 - lda [argp],Y - bpl cn0 - sec - lda #0 - sbc [argp] - sta [argp] - ldy #2 - lda #0 - sbc [argp],Y - sta [argp],Y - iny - iny - lda #0 - sbc [argp],Y - sta [argp],Y - iny - iny - lda #0 - sbc [argp],Y - sta [argp],Y - bra sn2 -sn0 lda ~isLong handle long values + lda ~isLong handle long and long long values beq sn0a ldy #2 - lda [argp],Y + lda ~isLongLong + beq sn0 + ldy #6 +sn0 lda [argp],Y bpl cn0 sec - lda #0 + ldx #0 + txa sbc [argp] sta [argp] - lda #0 + ldy #2 + txa + sbc [argp],Y + sta [argp],Y + lda ~isLongLong + beq sn2 + iny + iny + txa + sbc [argp],Y + sta [argp],Y + iny + iny + txa sbc [argp],Y sta [argp],Y bra sn2 @@ -3964,13 +3955,12 @@ cn1 lda [argp] push an int value cn1a pha cn2 ph4 #~str push the string addr ph2 #l:~str push the string buffer length - ph2 #0 do an unsigned conversion lda ~isLongLong do the proper conversion beq cn2a - pla jsr ~ULongLong2Dec bra pd1 -cn2a lda ~isLong +cn2a ph2 #0 do an unsigned conversion + lda ~isLong beq cn3 _Long2Dec bra pd1 @@ -4057,17 +4047,8 @@ pn1 lda ~hexPrefix if there is a hex prefix then jsl ~putchar ph2 ~hexPrefix+1 jsl ~putchar -pn1a lda ~paddChar if the number needs 0 padding then - cmp #'0' - bne pn1c - lda ~fieldWidth - bmi pn1c - beq pn1c -pn1b ph2 ~paddChar print padd zeros - jsl ~putchar - dec ~fieldWidth - bne pn1b -pn1c lda ~precision if the number needs more padding then +pn1a jsr ~ZeroPad pad with '0's if needed + lda ~precision if the number needs more padding then beq pn3 pn2 ph2 #'0' print padd characters jsl ~putchar @@ -4095,10 +4076,10 @@ pn5 cpy #l:~str quit if we're at the end of the ~str ; rn1 lda ~isLongLong beq rn2 - inc argp - inc argp - inc argp - inc argp + lda argp + clc + adc #4 + sta argp rn2 lda ~isLong beq rn3 inc argp @@ -4244,9 +4225,12 @@ lb1 clc restore the original argp+4 **************************************************************** * * ~Format_o - format an octal number +* ~Format_x - format a hexadecimal number (lowercase output) +* ~Format_X - format a hexadecimal number (uppercase output) +* ~Format_p - format a pointer * * Inputs: -* ~altForm - use a leading '0'? +* ~altForm - use a leading '0' (octal) or '0x' (hex)? * ~fieldWidth - output field width * ~paddChar - padd character * ~leftJustify - left justify the output? @@ -4260,15 +4244,34 @@ lb1 clc restore the original argp+4 ~Format_o private using ~printfCommon argp equ 7 argument pointer + + lda #3 use 3 bits per output character + bra cn0 + +~Format_x entry +; +; Set the "or" value; this is used to set the case of character results +; + lda #$20*256 + sta ~orVal + bra hx0 + +~Format_p entry + inc ~isLong +~Format_X entry + stz ~orVal +hx0 lda #4 use 4 bits per output character + ; ; Initialization ; +cn0 sta bitsPerChar + stz ~hexPrefix assume we won't lead with 0x stz ~sign ignore the sign flag lda #' ' initialize the string to blanks sta ~str move ~str,~str+1,#l:~str-1 - stz ~num+2 get the value to convert - lda ~isLongLong + lda ~isLongLong get the value to convert beq cn1 ldy #6 lda [argp],Y @@ -4279,7 +4282,7 @@ argp equ 7 argument pointer sta ~num+4 cn1 lda ~isLong beq cn2 -cn1a ldy #2 + ldy #2 lda [argp],Y sta ~num+2 cn2 lda [argp] @@ -4287,57 +4290,71 @@ cn2 lda [argp] beq cn2a and #$00FF cn2a sta ~num + ldx bitsPerChar if doing hex format then + cpx #3 + beq cn2b + ldx ~altForm if alt form has been selected then + beq cn2b + ora ~num+2 if value is not 0 then + ora ~num+4 + ora ~num+6 + beq cn2b + lda #'X0' set hex prefix to '0X' or '0x' + ora ~orVal + sta ~hexPrefix ; ; Convert the number to an ASCII string ; - short I,M - ldy #l:~str-1 set up the character index -cn3 lda ~num+7 quit if the number is zero - ora ~num+6 - ora ~num+5 - ora ~num+4 - ora ~num+3 - ora ~num+2 - ora ~num+1 - ora ~num - beq al1 - lda #0 roll off 3 bits - ldx #3 -cn4 lsr ~num+7 - ror ~num+6 - ror ~num+5 +cn2b ldy #l:~str-1 set up the character index +cn3 lda #' 0' roll off 4 bits + ldx bitsPerChar +cn4 lsr ~num+6 ror ~num+4 - ror ~num+3 ror ~num+2 - ror ~num+1 ror ~num ror A dex bne cn4 - lsr A form a character - lsr A - lsr A - lsr A - lsr A - ora #'0' + xba form a character + ldx bitsPerChar +cn4a asl A + dex + bne cn4a + cmp #('9'+1)*256+' ' if the character should be alpha then + blt cn5 + clc + adc #7*256 adjust it + ora ~orVal +cn5 dey sta ~str,Y save the character - dey - bra cn3 + lda ~num+6 loop if the number is not zero + ora ~num+4 + ora ~num+2 + ora ~num + bne cn3 ; -; If a leading zero is required, be sure we include one +; If a leading '0x' is required, be sure we include one ; -al1 cpy #l:~str-1 include a zero if no characters have - beq al2 been placed in the string - lda ~altForm branch if no leading zero is required + lda bitsPerChar if doing octal format then + cmp #3 + bne al3 + lda ~altForm if alt form has been selected then beq al3 -al2 lda #'0' - sta ~str,Y -al3 long I,M + lda ~precision make sure precision is non-zero + bne al2 + inc ~precision +al2 lda #'0 ' if the result is not ' 0' then + cmp ~str+l:~str-2 + beq al3 + sta ~str-1,Y include a zero in the string ; ; Piggy back off of ~Format_d for output ; - stz ~hexPrefix don't lead with 0x - brl ~Format_IntOut +al3 brl ~Format_IntOut +; +; Local data +; +bitsPerChar ds 2 bits per output character end **************************************************************** @@ -4357,36 +4374,36 @@ al3 long I,M using ~printfCommon argp equ 7 argument pointer - ph4