Geometry

[Andre Leao]

Released...



https://planetasinclair.blogspot.com/20 ... metry.html

✓ PROCESSED

DH 23/02/24. Added after the discussion and any updates that might have happened.

[ParadigmShifter]

The example screenshot I don't understand?

https://blogger.googleusercontent.com/i ... s320/2.jpg

What is perimeter1 and perimeter2? There's no simple formula for the perimeter of an ellipse... it's the very famous elliptic integral which not even Wolfram Alpha simplifies

(EDIT: See https://www.wolframalpha.com/input?i=el ... D34+r%3D12

also see

https://mathworld.wolfram.com/CompleteE ... dKind.html

Yikes)

https://www.mathsisfun.com/geometry/ell ... meter.html

It's approx 152.9 though which is neither of the results given there

[zarsoft]

Approximation means it is not correct.

My goal was to calculate an approximation with a simple formula
and not an approximation with good results.

[ParadigmShifter]

Ok well there's much more accurate approximation formulae on the mathsisfun site link.

I still don't understand what perimeter1 and perimeter2 are? Is the first one meant to be pi*R*r instead of pi*(R+r)? EDIT: No, pi*(R+r) makes more sense if it is nearly circular, my bad. So perimeter1 is the approximation if R and r are nearly the same (which they clearly aren't in this case)

It's kind of important to get maths stuff correct since there is only 1 "correct" answer to what is the perimeter of an ellipse (although it is hard to calculate). Both your approximations are way off though (first one because it isn't nearly circular so that's obvious). I think the second one is not good maybe needs to be 2*pi*SQRT((R^2+r^2)/2) ~= 160.19 that's approximation #1 on the mathsisfun link.

EDIT2: I think they have mistakenly tried to bring the /2 out of the sqrt in a way which is not correct.

EDIT3: Correct way to move the sqrt out the bracket is (pi/sqrt(2))*sqrt(R^2+r^2) or (1/sqrt(2))*pi*sqrt(R^2+r^2) = sqrt(1/2)*pi*sqrt(R^2+r^2). But the one with the /2 inside the brackets saves calculating sqrt(2) - but if you have sqrt(2) as a constant in the code use that instead. 1/sqrt(2) = sqrt(1/2) is probably another useful constant to have available.

You may want to watch this as well

[_dw]

The example screenshot I don't understand?

https://blogger.googleusercontent.com/i ... s320/2.jpg

What is perimeter1 and perimeter2? There's no simple formula for the perimeter of an ellipse... it's the very famous elliptic integral which not even Wolfram Alpha simplifies

(EDIT: See https://www.wolframalpha.com/input?i=el ... D34+r%3D12

also see

https://mathworld.wolfram.com/CompleteE ... dKind.html

Yikes)

https://www.mathsisfun.com/geometry/ell ... meter.html

It's approx 152.9 though which is neither of the results given there

I don't understand it at all, but that doesn't mean I won't talk about it.

At first glance, it is clear that these are supposed to be some approximate formulas based on the calculation of the circumference of a circle.

perimeter = 2 * PI * R = PI * (2*R)

adjusted to

approximate circumference of ellipse = PI * (R+r)

which probably arose from the consideration that the circumference will be smaller than 2 * PI * R and larger than 2 * PI * r. So I make an average.

approximate circumference of the ellipse = PI * ((R*R+r*r)**0.5)

The second pattern is very interesting because I make a circle like

approximate circumference of the ellipse = PI * hypotenuse calculated by the Pythagorean theorem

I assume that the reasoning was that the hypotenuse will always be greater than or equal to the longest length, but never greater than twice the length? And smaller or equal to the sum of the parts.

They can count.

Formulas from the Internet tell me

approximate circumference of the ellipse = PI * (R+r+(2*(R*R+r*r)**0.5) / 2 = PI * (12+34+51)/2 = 152.35

but elementary school is not enough to derive the formula..

[ParadigmShifter]

I'd watch the standup maths video first if I were you. (I've not watched it for a while though so not too familiar with what he says). I can probably help a bit though although elliptic integrals were way out of my league (although I do have a degree in maths - they are advanced stuff. They maybe did them in applied maths but I dropped that as soon as possible). Pure & stats only for me really

Degenerate case where R=1 r=0 has perimeter 4 (it's a line segment of length 2, which has "perimeter" = twice the length) which does not involve pi so there's a non-linear scaling of the factor of pi along the way as the differences between the major and minor radii change.

[_dw]

I think it's absolutely brilliant because it turns your brain on to deal with it. Because you get the feeling... I could prove this too!
And better!
I should be able to imagine it as an extremely elongated ellipse and count it as a circle with a smaller radius and take the longer part as segments. Make the ellipse a segment with rounded ends.

With the fact that the result will be the same or greater.

approximate circumference of ellipse = 2 * ( PI * r + 4 * (R - r)) = 2 * ( PI * 12 + 4 * 12 ) = 2 * ( PI * 12 + 12 ) = 171.4

Because when you unpack the integral on someone, I remember that in Saudi Arabia? there is a law that punishes witchcraft. One should be careful.

[ParadigmShifter]

It's not easy to prove.

It's basically calculation of an integral (which can be done numerically by computer to any desired accuracy). (length of a curve is an integral but they are usually difficult to solve analytically (analytically usually means in terms of elementary functions - "elementary" means that the result is defined in terms of elementary functions like polynomials, sin, cos and their inverses, exp, log, powers and roots), much harder than just doing a typical integral which calculates the area under a curve, not arc length).

Integrals (integration - area under a curve) is calculus, and integrals are much harder than derivatives (differentiation - slope of a curve at any point). Integrals hard whereas derivatives are easy usually.

Might want to start looking at this - if you know calculus (especially integration)

https://en.wikipedia.org/wiki/Arc_lengt ... ntegration

The calculation for a circle isn't that hard

https://math.stackexchange.com/question ... ntegration

but for an ellipse it is "way hard"

EDIT: I'm pretty sure the study of elliptic integrals (to calculate the arc length of ellipses) led to the study of "elliptic curves" - which don't look like ellipses at all btw, but they are involved with solving the elliptic integral in some way I don't understand - postgrad pure maths stuff. Elliptic curves were used by Andrew Wiles to prove Fermat's Last Theorem though That proof is extremely difficult to understand and is massive of course, proper PhD+++ level stuff.

EDIT2: Development of the elliptic integral

https://en.wikipedia.org/wiki/Elliptic_integral

[ParadigmShifter]

I think I missed out a factor of 2 in the equations I posted earlier for approx #1.

Math is hard - Barbie

2*(1/sqrt(2))*pi*sqrt(R^2+r^2) = (2/sqrt(2))*pi*sqrt(R^2+r^2) = 2*sqrt(1/2)*pi*sqrt(R^2+r^2) ~= 160.19 for R = 34, r = 12

I think that is the approximation you meant to use for perimeter2 anyway. (Which is 2pi * geometric mean of R and r, I think?). EDIT: Nope it's the root mean square of R and r

EDIT: That gives a value accurate to around 5% error if 3r < R

[Daveysloan]

This is exactly the sort of thing I'd show my parents to convince them to get a Spectrum, then spend all day playing Bruce Lee.

[Andre Leao]

It sounds so familiar

[zarsoft]

I think they have mistakenly tried to bring the /2 out of the sqrt in a way which is not correct.

You are right.

I've made that mistake when I was writing the program.

The correct formula should be
2*pi*sqrt((R^2+r^2)/2)

[ParadigmShifter]

Here's a circle drawing routine I wrote in ASM which you are welcome to use if you want

Code: Select all

	ORG #8000

RETURN_TO_BASIC EQU 1
NOPS_AT_START EQU 0
SCRBUF_BASEADDR EQU #4000
SPECIAL_CASE_LINE256 EQU 0

    MACRO ldim regpair, val1, val2
    ld regpair, ((val1&#FF)<<8)|(val2&#FF)
    ENDM

codestart:
    IF NOPS_AT_START
    nop
    nop
    nop
    nop
    ENDIF

main:
    IF RETURN_TO_BASIC
    ld (stashed_iy), iy
    exx
    ld (stashed_hl_alt), hl
	ld (stashed_sp), sp
    ENDIF

.mainloop
	ld a, 7
	out (#fe), a	
	halt
	ld a, 1
	out (#fe), a	

	IF 1
	ldim hl, 128, 96
	ld d, 63
	call draw_circle
	jp .mainloop
	ENDIF
	IF 0;1

	ldim hl, 128, 96
	ld d, 95
	call draw_circle

	ld b, 50
.morehalt
	halt
	djnz .morehalt
	ld a, 0|(7<<3)
	call cls

	ldim hl, 128, 96
	ld d, 95
	call draw_fill_circle

	ld b, 50
.morehalt2
	halt
	djnz .morehalt2
	ld a, 0|(7<<3)
	call cls

	jp .mainloop
	ENDIF
	IF 1;0
	ldim hl, 128, 96
	ld d, 95
.nextrad
	push de
	push hl
	call draw_circle
	pop hl
	pop de
	dec d
	jp p, .nextrad
	ENDIF

    IF RETURN_TO_BASIC
returntobasic:
    ld iy, (stashed_iy)
    ld hl, (stashed_hl_alt)
	ld sp, (stashed_sp)
    exx
    ENDIF
	ret

; H - cx
; L - cy
; D - radius
draw_circle:
	; B = xoffset
	xor a
	; A' = sx
	ex af, af'
	xor a
	; C = yoffset
	ld b, a
	ld c, d
	ld d, a

	IF 1
	; push everything
	; do the very first line, just draw 2 pixels
	push de
	push bc
	push hl
	; cy+x
	;ld a, l
	;add b
	;ld l, a ; B is 0 first line

	; cx-y
	ld a, h
	sub c
	ld b, a

	; cx+y
	ld a, h
	add c
	ld c, a

	xor a ; A=0 and clear carry
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l
	ld h, (hl)
	ld l, a

	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, c
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a

	ld a, l
	and ~31
	ld l, a
	ld a, b
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, b
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a
	jp .doneveryfirstline
	ENDIF
.nextline

	push de
	push bc
	push hl
	; cy+x
	ld a, l
	add b
	ld l, a

	; cx-y
	ld a, h
	sub c
	ld b, a

	; cx+y
	ld a, h
	add c
	ld c, a

	xor a ; A=0 and clear carry
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l
	ld h, (hl)
	ld l, a

	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, c
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a

	ld a, l
	and ~31
	ld l, a
	ld a, b
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, b
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a

	pop hl
	pop bc
	push bc
	push hl
	; cy-x
	ld a, l
	sub b
	ld l, a

	; cx-y
	ld a, h
	sub c
	ld b, a

	; cx+y
	ld a, h
	add c
	ld c, a

	xor a ; A=0 and clear carry
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l
	ld h, (hl)
	ld l, a

	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, c
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a

	ld a, l
	and ~31
	ld l, a
	ld a, b
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, b
	ld de, col2pix
	and 7
	add e
	ld e, a
	ld a, (de)
	or (hl)
	ld (hl), a

.doneveryfirstline
	pop hl
	pop bc
	pop de

	push hl
	ld h, 0
	ld l, d
	ld d, h
	ld e, b
	add hl, de
	add hl, de
	inc hl
	inc b
	ld a, l ; A is new D if we jump to .nexty soon
	ld d, 0
	ld e, c
	or a
	sbc hl, de
	or a
	sbc hl, de
	inc hl
	bit 7, h
	ld d, a
	jr nz, .nexty
	ld d, l

	pop hl
	ld a, b
	sub c
	dec a
	jr z, .justdecy ; we already drew this line (at 45 degrees)

	push de

	push bc
	push hl
	; cy-y
	ld a, l
	sub c
	ld l, a
	; cx - x + 1
	ex af, af'
	ld e, a
	ex af, af'
	ld a, h
	sub b
	inc a
	ld c, a
	; b = sx - x
	ld a, b
	sub e
	ld b, a

	xor a
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l ; this won't overflow because of align 256
	ld h, (hl)
	ld l, a

	push hl
	call draw_line_horz_knowaddr
	pop de

	pop hl
	pop bc
	push bc
	push hl
	; cx + sx
	ld a, h
	ex af, af'
	ld h, a
	ex af, af'
	add h
	ld c, a
	; b = x - sx
	ld a, b
	sub h
	ld b, a

	ld h, d
	ld l, e
	call draw_line_horz_knowaddr
	pop hl
	pop bc

	push bc
	push hl
	; cy+y
	ld a, l
	add c
	ld l, a
	; cx - x + 1
	ld a, h
	sub b
	inc a
	ld c, a
	; b = sx - x
	ld a, b
	ex af, af'
	ld b, a
	ex af, af'
	sub b
	ld b, a

	xor a
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l ; this won't overflow because of align 256
	ld h, (hl)
	ld l, a

	push hl
	call draw_line_horz_knowaddr
	pop de
	pop hl
	pop bc

	push bc
	push hl
	; cx + sx
	ld a, h
	ex af, af'
	ld h, a
	ex af, af'
	add h
	ld c, a
	; b = x - sx
	ld a, b
	sub h
	ld b, a

	ld h, d
	ld l, e
	call draw_line_horz_knowaddr
	pop hl

	pop bc

	pop de

.justdecy
	dec c

	ex af, af'
	ld a, b
	ex af, af'

	push hl
.nexty
	ld a, c
	sub b
	pop hl
	jp p, .nextline
	ret

; H - cx
; L - cy
; D - radius
draw_fill_circle:
	; B = xoffset
	xor a
	; C = yoffset
	ld b, a
	ld c, d
	ld d, a
	; push everything
	push de

	; draw the very first line
	push bc
	push hl
	; cy+x
	;ld a, l
	;add b
	;ld l, a ; B is 0 here
	ld a, c
	add c
	inc a
	ld b, a
	; cx-y
	ld a, h
	sub c
	ld c, a
	call draw_line_horz_xor
	jp .doneveryfirstline
.nextline
	push de

	push bc
	push hl
	; cy+x
	ld a, l
	add b
	ld l, a
	ld a, c
	add c
	inc a
	ld b, a
	; cx-y
	ld a, h
	sub c
	ld c, a
	call draw_line_horz_xor
	pop hl
	pop bc
	push bc
	push hl
	; cy-x
	ld a, l
	sub b
	ld l, a
	ld a, c
	add c
	inc a
	ld b, a
	; cx-y
	ld a, h
	sub c
	ld c, a
	call draw_line_horz_xor

.doneveryfirstline
	pop hl
	pop bc

	; work out new d
	pop de

	push hl
	ld h, 0
	ld l, d
	ld d, h
	ld e, b
	add hl, de
	add hl, de
	inc hl
	inc b
	ld a, l ; A is new D if we jump to .nexty soon
	ld d, 0
	ld e, c
	or a
	sbc hl, de
	or a
	sbc hl, de
	inc hl
	bit 7, h
	ld d, a
	jr nz, .nexty

	ld d, l
	pop hl
	ld a, b
	sub c
	dec a
	jr z, .justdecy ; we already drew this line (at 45 degrees)

	push de

	push bc
	push hl
	; cy-y
	ld a, l
	sub c
	ld l, a
	; cx - x + 1
	ld a, h
	sub b
	inc a
	ld c, a
	; b = x+x-1
	ld a, b
	add b
	dec a
	ld b, a
	call draw_line_horz_xor
	pop hl
	pop bc
	push bc
	push hl
	; cy+y
	ld a, l
	add c
	ld l, a
	; cx - x + 1
	ld a, h
	sub b
	inc a
	ld c, a
	; b = x+x-1
	ld a, b
	add b
	dec a
	ld b, a
	call draw_line_horz_xor
	pop hl
	pop bc

	pop de

.justdecy
	dec c

	push hl
.nexty
	ld a, c
	sub b
	pop hl
	jp p, .nextline
	ret

	IF 0
; C - column
; L - row
plot:
	xor a
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l
	ld h, (hl)
	ld l, a
	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	ld a, c
	exx
	ld hl, col2pix
	and 7
	add l
	ld l, a
	ld a, (hl)
	exx
	or (hl)
	ld (hl), a
	ret
	ENDIF

; B - number of pixels
; C - column offset
; L - row
draw_line_horz:
	xor a
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l ; this won't overflow because of align 256
	ld h, (hl)
	; work out the column address
	ld l, a
draw_line_horz_knowaddr:
	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	; if C&7 == 0 and B&7 == 0, just draw the middle of the line
	ld a, c
	or b
	and 7
	jr z, .drawlinemiddle
.notmiddleonly
	; work out (C&7), amount of pixels to draw on left hand side
	ld a, c
	and 7
	jr z, .justdrawrhs
	ld e, a ; remember (C&7) in E. will subtract this from line length in a bit
	; if 8-C&7 > B (line length), we only draw 1 cell and not all of the pixels
	ld a, 8
	sub e
	cp b
	jr nc, .draw1cellpartial
	ld c, a
	ld a, e ; amount to shift down
	exx
	ld hl, srlFFtable
	add l
	ld l, a
	ld a, (hl)
	exx
	; put lefthandside
	or (hl)
	ld (hl), a
	inc l
	ld a, b
	sub c ; adjust line length for pixels we just drew
	ld b, a
	and 7
.justdrawrhs
	ld e, b ; remember B&7 in E
	ld a, b
	and ~7
	rrca
	rrca
	rrca
	ld b, a
	or a
	jr z, .dontdrawmiddle
	ld a, #FF
.middleloop
	ld (hl), a
	inc l
	djnz .middleloop
.dontdrawmiddle
	; work out (B&7), amount of pixels to draw on left hand side
	ld a, e
	and 7
	dec a
	ret m
	jr z, .putlinerhsnoshift ; no need to shift
	exx
	ld hl, sra80table
	add l
	ld l, a
	ld a, (hl)
	exx
	; put righthandside
	or (hl)
	ld (hl), a
	ret
.putlinerhsnoshift
	ld a, #80
	or (hl)
	ld (hl), a
	ret
.draw1cellpartial
	dec b
	ld a, #80
	jr z, .nosra
	ld a, b
	exx
	ld hl, sra80table
	add l
	ld l, a
	ld a, (hl)
	exx
.nosra
	ld b, e
.rrcaagain
	rrca
	djnz .rrcaagain
	or (hl)
	ld (hl), a
	ret
.drawlinemiddle
	ld a, b
	rrca
	rrca
	rrca
	ld b, a
	ld a, #FF
.middleonlyloop
	ld (hl), a
	inc l
	djnz .middleonlyloop
	ret

; B - number of pixels
; C - column offset
; L - row
draw_line_horz_xor:
	xor a
	sla l
	ld h, tbl_scraddr/256 ; this requires tbl_scraddr be aligned to a 256 byte boundary
	adc h ; A = H + carry flag
	ld h, a ; write it back
	; look up screen address from  the table
	ld a, (hl)
	inc l ; this won't overflow because of align 256
	ld h, (hl)
	; work out the column address
	ld l, a
	ld a, c
	and ~7
	rrca
	rrca
	rrca
	add l
	ld l, a
	; if C&7 == 0 and B&7 == 0, just draw the middle of the line
	ld a, c
	or b
	and 7
	jr z, .drawlinemiddle
.notmiddleonly
	; work out (C&7), amount of pixels to draw on left hand side
	ld a, c
	and 7
	jr z, .justdrawrhs
	ld e, a ; remember (C&7) in E. will subtract this from line length in a bit
	; if 8-C&7 > B (line length), we only draw 1 cell and not all of the pixels
	ld a, 8
	sub e
	cp b
	jr nc, .draw1cellpartial
	ld c, a
	ld a, e ; amount to shift down
	exx
	ld hl, srlFFtable
	add l
	ld l, a
	ld a, (hl)
	exx
	; put lefthandside
	xor (hl)
	ld (hl), a
	inc l
	ld a, b
	sub c ; adjust line length for pixels we just drew
	ld b, a
	and 7
.justdrawrhs
	ld e, b ; remember B&7 in E
	ld a, b
	and ~7
	rrca
	rrca
	rrca
	ld b, a
	or a
	jr z, .dontdrawmiddle
	;ld a, #FF
	ld c, #FF
.middleloop
	ld a, c
	xor (hl)
	ld (hl), a
	inc l
	djnz .middleloop
.dontdrawmiddle
	; work out (B&7), amount of pixels to draw on left hand side
	ld a, e
	and 7
	dec a
	ret m
	jr z, .putlinerhsnoshift ; no need to shift
	exx
	ld hl, sra80table
	add l
	ld l, a
	ld a, (hl)
	exx
	; put righthandside
	xor (hl)
	ld (hl), a
	ret
.putlinerhsnoshift
	ld a, #80
	xor (hl)
	ld (hl), a
	ret
.draw1cellpartial
	dec b
	ld a, #80
	jr z, .nosra
	ld a, b
	exx
	ld hl, sra80table
	add l
	ld l, a
	ld a, (hl)
	exx
.nosra
	ld b, e
.rrcaagain
	rrca
	djnz .rrcaagain
	xor (hl)
	ld (hl), a
	ret
.drawlinemiddle
	ld a, b
	rrca
	rrca
	rrca
	ld b, a
	ld c, #FF
.middleonlyloop
	ld a, c
	xor (hl)
	ld (hl), a
	inc l
	djnz .middleonlyloop
	ret

cls:
	di                  ;disable interrupt
	ld (.stack+1), sp	;store current stack pointer
	ld sp, 16384 + 6144 + 768
	ld b, 128			; clear attribs in 128 * 3 pushes
	ld h, a
	ld l, a
.attribloop
	push hl
	push hl
	push hl
	djnz .attribloop
	ld hl, 0
	ld b, l             ;set B to 0. it causes that DJNZ will repeat 256 times
.loop1
	push hl             ;store hl on stack
	push hl             ;next
	push hl             ;these four push instruction stores 8 bytes on stack
	push hl
	push hl             ;store hl on stack
	push hl             ;next
	push hl             ;these four push instruction stores 8 bytes on stack
	push hl
	push hl             ;store hl on stack
	push hl             ;next
	push hl             ;these four push instruction stores 8 bytes on stack
	push hl
	djnz .loop1			;repeat for next 12*2 bytes
.stack
	ld sp, 0            ;parameter will be overwritten
	ei
	ret

data_section:
    IF RETURN_TO_BASIC
stashed_iy dw 0
stashed_hl_alt dw 0
stashed_sp dw 0
    ENDIF

	ALIGN 256
; screen address table. This must be 256 byte aligned
tbl_scraddr	dw SCRBUF_BASEADDR + #0000, SCRBUF_BASEADDR + #0100, SCRBUF_BASEADDR + #0200, SCRBUF_BASEADDR + #0300, SCRBUF_BASEADDR + #0400, SCRBUF_BASEADDR + #0500, SCRBUF_BASEADDR + #0600, SCRBUF_BASEADDR + #0700 
			dw SCRBUF_BASEADDR + #0020, SCRBUF_BASEADDR + #0120, SCRBUF_BASEADDR + #0220, SCRBUF_BASEADDR + #0320, SCRBUF_BASEADDR + #0420, SCRBUF_BASEADDR + #0520, SCRBUF_BASEADDR + #0620, SCRBUF_BASEADDR + #0720
			dw SCRBUF_BASEADDR + #0040, SCRBUF_BASEADDR + #0140, SCRBUF_BASEADDR + #0240, SCRBUF_BASEADDR + #0340, SCRBUF_BASEADDR + #0440, SCRBUF_BASEADDR + #0540, SCRBUF_BASEADDR + #0640, SCRBUF_BASEADDR + #0740
			dw SCRBUF_BASEADDR + #0060, SCRBUF_BASEADDR + #0160, SCRBUF_BASEADDR + #0260, SCRBUF_BASEADDR + #0360, SCRBUF_BASEADDR + #0460, SCRBUF_BASEADDR + #0560, SCRBUF_BASEADDR + #0660, SCRBUF_BASEADDR + #0760
			dw SCRBUF_BASEADDR + #0080, SCRBUF_BASEADDR + #0180, SCRBUF_BASEADDR + #0280, SCRBUF_BASEADDR + #0380, SCRBUF_BASEADDR + #0480, SCRBUF_BASEADDR + #0580, SCRBUF_BASEADDR + #0680, SCRBUF_BASEADDR + #0780
			dw SCRBUF_BASEADDR + #00a0, SCRBUF_BASEADDR + #01a0, SCRBUF_BASEADDR + #02a0, SCRBUF_BASEADDR + #03a0, SCRBUF_BASEADDR + #04a0, SCRBUF_BASEADDR + #05a0, SCRBUF_BASEADDR + #06a0, SCRBUF_BASEADDR + #07a0
			dw SCRBUF_BASEADDR + #00c0, SCRBUF_BASEADDR + #01c0, SCRBUF_BASEADDR + #02c0, SCRBUF_BASEADDR + #03c0, SCRBUF_BASEADDR + #04c0, SCRBUF_BASEADDR + #05c0, SCRBUF_BASEADDR + #06c0, SCRBUF_BASEADDR + #07c0
			dw SCRBUF_BASEADDR + #00e0, SCRBUF_BASEADDR + #01e0, SCRBUF_BASEADDR + #02e0, SCRBUF_BASEADDR + #03e0, SCRBUF_BASEADDR + #04e0, SCRBUF_BASEADDR + #05e0, SCRBUF_BASEADDR + #06e0, SCRBUF_BASEADDR + #07e0
			dw SCRBUF_BASEADDR + #0800, SCRBUF_BASEADDR + #0900, SCRBUF_BASEADDR + #0a00, SCRBUF_BASEADDR + #0b00, SCRBUF_BASEADDR + #0c00, SCRBUF_BASEADDR + #0d00, SCRBUF_BASEADDR + #0e00, SCRBUF_BASEADDR + #0f00 
			dw SCRBUF_BASEADDR + #0820, SCRBUF_BASEADDR + #0920, SCRBUF_BASEADDR + #0a20, SCRBUF_BASEADDR + #0b20, SCRBUF_BASEADDR + #0c20, SCRBUF_BASEADDR + #0d20, SCRBUF_BASEADDR + #0e20, SCRBUF_BASEADDR + #0f20
			dw SCRBUF_BASEADDR + #0840, SCRBUF_BASEADDR + #0940, SCRBUF_BASEADDR + #0a40, SCRBUF_BASEADDR + #0b40, SCRBUF_BASEADDR + #0c40, SCRBUF_BASEADDR + #0d40, SCRBUF_BASEADDR + #0e40, SCRBUF_BASEADDR + #0f40
			dw SCRBUF_BASEADDR + #0860, SCRBUF_BASEADDR + #0960, SCRBUF_BASEADDR + #0a60, SCRBUF_BASEADDR + #0b60, SCRBUF_BASEADDR + #0c60, SCRBUF_BASEADDR + #0d60, SCRBUF_BASEADDR + #0e60, SCRBUF_BASEADDR + #0f60
			dw SCRBUF_BASEADDR + #0880, SCRBUF_BASEADDR + #0980, SCRBUF_BASEADDR + #0a80, SCRBUF_BASEADDR + #0b80, SCRBUF_BASEADDR + #0c80, SCRBUF_BASEADDR + #0d80, SCRBUF_BASEADDR + #0e80, SCRBUF_BASEADDR + #0f80
			dw SCRBUF_BASEADDR + #08a0, SCRBUF_BASEADDR + #09a0, SCRBUF_BASEADDR + #0aa0, SCRBUF_BASEADDR + #0ba0, SCRBUF_BASEADDR + #0ca0, SCRBUF_BASEADDR + #0da0, SCRBUF_BASEADDR + #0ea0, SCRBUF_BASEADDR + #0fa0
			dw SCRBUF_BASEADDR + #08c0, SCRBUF_BASEADDR + #09c0, SCRBUF_BASEADDR + #0ac0, SCRBUF_BASEADDR + #0bc0, SCRBUF_BASEADDR + #0cc0, SCRBUF_BASEADDR + #0dc0, SCRBUF_BASEADDR + #0ec0, SCRBUF_BASEADDR + #0fc0
			dw SCRBUF_BASEADDR + #08e0, SCRBUF_BASEADDR + #09e0, SCRBUF_BASEADDR + #0ae0, SCRBUF_BASEADDR + #0be0, SCRBUF_BASEADDR + #0ce0, SCRBUF_BASEADDR + #0de0, SCRBUF_BASEADDR + #0ee0, SCRBUF_BASEADDR + #0fe0
			dw SCRBUF_BASEADDR + #1000, SCRBUF_BASEADDR + #1100, SCRBUF_BASEADDR + #1200, SCRBUF_BASEADDR + #1300, SCRBUF_BASEADDR + #1400, SCRBUF_BASEADDR + #1500, SCRBUF_BASEADDR + #1600, SCRBUF_BASEADDR + #1700
			dw SCRBUF_BASEADDR + #1020, SCRBUF_BASEADDR + #1120, SCRBUF_BASEADDR + #1220, SCRBUF_BASEADDR + #1320, SCRBUF_BASEADDR + #1420, SCRBUF_BASEADDR + #1520, SCRBUF_BASEADDR + #1620, SCRBUF_BASEADDR + #1720
			dw SCRBUF_BASEADDR + #1040, SCRBUF_BASEADDR + #1140, SCRBUF_BASEADDR + #1240, SCRBUF_BASEADDR + #1340, SCRBUF_BASEADDR + #1440, SCRBUF_BASEADDR + #1540, SCRBUF_BASEADDR + #1640, SCRBUF_BASEADDR + #1740
			dw SCRBUF_BASEADDR + #1060, SCRBUF_BASEADDR + #1160, SCRBUF_BASEADDR + #1260, SCRBUF_BASEADDR + #1360, SCRBUF_BASEADDR + #1460, SCRBUF_BASEADDR + #1560, SCRBUF_BASEADDR + #1660, SCRBUF_BASEADDR + #1760
			dw SCRBUF_BASEADDR + #1080, SCRBUF_BASEADDR + #1180, SCRBUF_BASEADDR + #1280, SCRBUF_BASEADDR + #1380, SCRBUF_BASEADDR + #1480, SCRBUF_BASEADDR + #1580, SCRBUF_BASEADDR + #1680, SCRBUF_BASEADDR + #1780
			dw SCRBUF_BASEADDR + #10a0, SCRBUF_BASEADDR + #11a0, SCRBUF_BASEADDR + #12a0, SCRBUF_BASEADDR + #13a0, SCRBUF_BASEADDR + #14a0, SCRBUF_BASEADDR + #15a0, SCRBUF_BASEADDR + #16a0, SCRBUF_BASEADDR + #17a0
			dw SCRBUF_BASEADDR + #10c0, SCRBUF_BASEADDR + #11c0, SCRBUF_BASEADDR + #12c0, SCRBUF_BASEADDR + #13c0, SCRBUF_BASEADDR + #14c0, SCRBUF_BASEADDR + #15c0, SCRBUF_BASEADDR + #16c0, SCRBUF_BASEADDR + #17c0
			dw SCRBUF_BASEADDR + #10e0, SCRBUF_BASEADDR + #11e0, SCRBUF_BASEADDR + #12e0, SCRBUF_BASEADDR + #13e0, SCRBUF_BASEADDR + #14e0, SCRBUF_BASEADDR + #15e0, SCRBUF_BASEADDR + #16e0, SCRBUF_BASEADDR + #17e0

col2pix	db #80, #40, #20, #10, #08, #04, #02, #01
; #FF >> 0..7
srlFFtable db #FF, #7F, #3F, #1F, #0F, #07, #03, #01
; same but for arithmetic shift right
sra80table db #80, #C0, #E0, #F0, #F8, #FC, #FE, #FF
endofprog:

upload jpg

It only draws circles on screen... I got y-clipping going but x-clipping was a bit gnarly.

It's based on this basic code which can easily be adapted to draw an ellipse. It's an implementation of this algorithm

https://en.wikipedia.org/wiki/Midpoint_circle_algorithm

which is widely regarded as the best circle drawing algorithm since it only uses addition and subtraction (and you can also optimise x+x using a shift which I did in the ASM).

I think this draws filled circles though... the code to draw unfilled circles is commented out

Code: Select all

5 REM cx,cy are the circle centre, r is the radius
10 LET cx=128: LET cy=88: LET r=80
15 REM x, y are the offsets from the centre
16 REM d is the difference between the distance from the centre
17 REM squared, and the actual radius squared
20 LET x=0: LET y=r: LET d=0
25 REM plot all eight quadrants
30 REM PLOT cx+x,cy+y
40 REM PLOT cx+x,cy-y
50 PLOT cx-x,cy+y: DRAW x+x,0
60 PLOT cx-x,cy-y: DRAW x+x,0
70 REM PLOT cx+y,cy+x
80 PLOT cx-y,cy+x: DRAW y+y,0
90 REM PLOT cx+y,cy-x
100 PLOT cx-y,cy-x: DRAW y+y,0
105 REM compute the new d if we increase x
106 REM d=(x*x)+(y*y)-(r*r)
107 REM d2=(x+1)*(x+1)+(y*y)-(r*r)=(x*x)+(2*x)+1+(y*y)-(r*r)
108 REM d2-d = (2*x)+1
110 LET d=d+x+x+1
120 LET x=x+1
125 REM check if we should move closer to the circle
126 REM e=new d if y=y-1 using similar equation as above
130 LET e=d-y-y+1
140 IF e<0 THEN GO TO 170
150 LET d=e
160 LET y=y-1
165 REM once y<x we are past 45 degrees, so stop
170 IF y>=x THEN GO TO 30
175 REM keep drawing smaller circles down to zero
180 STOP : IF r=0 THEN STOP 
190 LET r=r-1
200 GO TO 20

I also optimised the non-filled version to draw horizontal line segments instead of individual pixels (which the machine code version already does), but the BASIC code was posted in a thread on WoS which was then locked and hidden so that is lost now, wasn't too hard to do though.

[_dw]

You are right.

I've made that mistake when I was writing the program.

The correct formula should be
2*pi*sqrt((R^2+r^2)/2)

I would think it over once more because during the calculation
PI * (R+r) is felt from the diameter of the circumference of the circle of the calculated larger number and the circumference of the circle of the smaller number. And it's a nice formula for R = approximately r. Both circuits form an "envelope" around the actual result.

But if you calculate it from the hypotenuse, it will always be greater (or equal when r=0) than R. So greater than R, which forms the outer/larger "envelope".

[ParadigmShifter]

That formula is the correct (intended) one (it's approximation #1 on mathsisfun link) and is the first decent one Matt Parker mentions in the video I posted. I think you are right that it always overestimates the perimeter though.

Dunno what you mean by the "envelope", sorry... Maybe you mean the convex hull or something.

pi*(R+r) is terrible if R and r are not of similar magnitude...

The sqrt((R^2+r^2)/2) is the root mean square error of R and r

https://en.wikipedia.org/wiki/Root_mean_square

A much much better formula is the one from Ramanujan mentioned in Matt's video

C ~= pi*(3*(R+r) - sqrt((3*R+r)(R+3*r)))

and the one he mentions next which uses the value h = ((R-r)^2)/((R+r)^2) is even better

C ~= pi*(R+r)*(1+(3*h/(10+sqrt(4-3*h))))

see https://en.wikipedia.org/wiki/Ellipse#Circumference

Ramanujan was of course a genius

[zarsoft]

I would think it over once...
Both circuits form an "envelope" around the actual result.

OK, then I will use P3 = (P1+P2)/2

Approximation 1:
P1 = pi*(R+r)
P1 = pi*(34+12) = 144.5

Approximation 2:
P2 = 2*pi*sqrt((R^2+r^2)/2)
P2 = 2*pi*sqrt((34^2+12^2)/2) = 160.19

Approximation 3:
P3 = (P1+P2)/2
P3 = (144.5+160.19) = 152.345

This is a reasonably approximation.

Of course there are better options...

Approximation 4:
P4 = pi*(3*(R+r)−sqrt((3*R+r)*(R+3*r)))
P4 = pi*(3*(34+12)−sqrt((3*34+12)*(34+3*12))) = 152.8987

Approximation 5:
(Best approximation, widely used by engineers for serious work at the end of the 20th century)
L = (R-r)/(R+r)
P5 = pi*(R+r)*(1+(3*L^2)/(10+sqrt(4-3*L^2)))
P5 = pi*(34+12)*(1+(3*((34-12)/(34+12))^2)/(10+sqrt(4-3*((34-12)/(34+12))^2))) = 152.9026

[_dw]

By "envelope" I meant the term where you get a guaranteed value greater/lesser (in extreme cases or the same) than the one you are looking for.

Now I see that it is not a hypotenuse because it is already divided by 2. I missed that.

[ParadigmShifter]

Those last 2 are Ramanujan's approximations. You're just using sqrt(h) and then squaring in #5.

Those are the "best" simple ones which don't involve terms of a power series.

One of the power series ones converges very quickly and all the denominators are powers of 2 which is nice.

"Infinite Series 2" from here

https://www.mathsisfun.com/geometry/ell ... meter.html

which uses falling factorials

https://en.wikipedia.org/wiki/Binomial_ ... ial_series

https://en.wikipedia.org/wiki/Falling_a ... factorials

watch out for the "alpha to the k underlined" notation which is the falling factorial



So the terms in the power series aren't hard to work out and you can put them in a lookup table.

[ParadigmShifter]

By "envelope" I meant the term where you get a guaranteed value greater/lesser (in extreme cases or the same) than the one you are looking for.

Now I see that it is not a hypotenuse because it is already divided by 2. I missed that.

Ok you mean the error term. Root mean square is a popular way of estimating the error and works pretty well in this case (<5% error if 3r<R).


EDIT: It's probably worth stating that knowing the exact value (or a good approximation) of the perimeter of an ellipse is not that useful unless you want to work out the lengths of planetary orbits I suppose.

Ellipses and ellipsoids aren't used much in computer games (so for collision detection and physics) because working out ellipse/ellipse intersections is very difficult (involves solving a quartic polynomial).

The equation for solving quartics is very long...



... which is why most games use simpler collision primitives like cylinders and capsules if they want something approximately ellipse/ellipsoid shaped.

However, doing ellipse vs triangle mesh collision is pretty easy (you just scale the ellipse back up to a circle or sphere and then scale the triangle mesh accordingly).

[Andre Leao]

Zé Oliveira fixed the program. You can download the new version here:

https://planetasinclair.blogspot.com/20 ... igido.html

[ParadigmShifter]

Using the average of the 2 approximations (although the circle approximation pi*R*r is really not good) does seem to improve the result of the 2nd approximation so I am wondering why that was not mentioned in any of the places I have seen.

I am guessing that is because it seems to always produce an underestimate from some back of the envelope tests so isn't as useful in real world situations (probably engineering, what is the length of kerbstones required to fit round an elliptical roundabout) so it's generally better to have some material left over than not enough? And Ramanujan's approximations are probably better anyway? (another guess there though - Ramanujan was a brilliant genius though so I would have thought he would have considered the average method himself at first?)

I'd suggest just using the first Ramanujan approximation (the sqrt((3R+r)(R+3r)) one) myself.

I'd also introduce the approx equal sign as well (what I am using ~= for but should actually be ≅ or ≈) since stating perimeter = (A1+ A2) / 2 is not correct really. All the approximations should use that sign instead of = since equality has a specific meaning of course and approximations should not use = sign It's ok to use approx1 = (approximation formula) and approx2 = (approximation formula) since those are the definitions of approx1 and approx2 but saying perimeter = (A1+A2)/2 is NOT correct.

Gonna have to read up a bit more on this I guess Or it might be worth asking the maths dudes who probably know more about the subject than myself on Matt Parker's video (since he was asking for better approximations) - although it is quite an old video I think so maybe won't be seen if we do that.

[ParadigmShifter]

Ok so this site

https://www.cuemath.com/measurement/per ... f-ellipse/

suggests averaging these 2 approximations

P ≈ π √[ 2 (a^2 + b^2) ]
P ≈ π [ (3/2)(a+b) - √(ab) ]

The first one is just your approximation2, (approximation1 on mathsisfun site) rearranged slightly

(see: https://www.wolframalpha.com/input?i=2* ... %2F2%29%29 and using the fact that 2/√2 = √2 (since that's (2^1)/(2^(1/2)) = 2^(1/2) = √2 using the identity: a^n/a^m = a^(n-m)), and then moving the √2 inside the square root using √a√b = √(ab))

P ≈ π √[ 2 (a^2 + b^2) ] is a nicer looking equation than we had before I think since it gets rid of the division.

I've not seen the second formula before though but it does mention it underestimates.



But again they say Ramanujan's approximations are better.

[zarsoft]

perimeter = (A1+A2)/2 is NOT correct.

P ≈ π √[ 2 (a^2 + b^2) ] is a nicer looking equation

You are right.
But now I'm not going to change the program anymore...

[ParadigmShifter]

No probs. It's good that you removed the algebra error though moving the /2 out of the brackets incorrectly. (I also made an error by forgetting about the factor of 2 at my first attempt of simplification as well and I failed to spot that 2/√2 = √2, which I should have realised since in the days when sqrt was slow on computers we used to use a fast reciprocal square root rsqrt(x) ~= 1/sqrt(x) approximation and have sqrt(x) return x*rsqrt(x). rsqrt is sometimes also called isqrt for inverse square root but I don't like that name since the inverse of square root is squaring

Reciprocal sqrt approximation used in Quake 3



Nowadays that "optimisation" is not good since sqrt(x) is more efficient on modern FPUs. It was good back in the day though (will not work on spectrum since it relies on IEEE floating point format which the speccy does not use).

https://en.wikipedia.org/wiki/Fast_inve ... solescence

[ParadigmShifter]

Ok I also won't moan about the PI symbol you use looks like either a capital pi (which would be the symbol for product in same way capital sigma is used for sum) or else suggests that pi is a vector :p

Next step for geometry would of course be a conversion of this excellent ruler and compass construction tool to the speccy

https://www.mathspad.co.uk/i2/construct.php

Although I studied those in Fields & Galois theory (ruler and compass constructions are a field extension of Q, the field of rational numbers) and you can use Galois Theory (which I was terrible at since it's very hard to understand) to show that trisecting an arbitrary angle, doubling the cube and squaring the circle are impossible to perform using classical ruler and compass constructions only. You can do some of those using a marked ruler though). EDIT: Looks like you don't need Galois Theory to prove that, just the theory of Field Extensions. It was a long time ago and I skived most of the lectures lol. Reason you can't do those things is because ruler and compass constructions only allow you to solve linear and quadratic equations (when using a collapsible compass i.e. you can't use it to remember a length and the ruler has no markings, so just a straight edge) whereas trisecting an arbitrary angle and doubling the cube both require constructing a cube root and squaring the circle would require constructing a length equal to the transcendental number pi). pi was not proven to be transcendental until 1882

I can still remember how to construct an equilateral triangle though



EDIT2: And of course you can construct a perpendicular bisector of a line segment by constructing another equilateral triangle on the other side of the line and connecting the points of both triangles not on the original line segment. Most youtube videos showing how to do that assume the compass does not collapse when lifted off the page which isn't allowed under classical Greek rules