| | Scout
Registered: Dec 2002
Posts: 1568 |
3D projection on the C=64 (or...how do I divide?)
Hi!
After 12 years I picked up coding again on the C=64.
What I want to make is a simple 3D starfield.
I know something about 3D programming, so that's no issue.
The only thing I stumbled on is how to implement 3D projection (X1=X/Z) and especially the dividing part of it.
How can I implement this on the C=64?
Any hints, tips?
Thanx!
R.
(Asm on the PC is so much easier ;-) |
|
| | Perff
Administrator
Posts: 1665 |
One of the methods I have used the most when dividing (and multiplying) in assembler is using logarithm.
I use the fact that:
x/y = e^(log x - log y)
You have one table containing values of log(x), and another table of e^x
When you want to divide the content of the X-reg with the Y-reg you do:
sec
lda logtable,x
sbc logtable,y
tax
lda exptable,x
That's it! :)
(You must of course take care of overflows etc.)
One downside using this is that the result is not 100% accurate, but it's fast and you can also use it for multiplying. (do an add instead of sbc) |
| | Scout
Registered: Dec 2002
Posts: 1568 |
Hey!
Now that's what I call a great solution!
Cheers!
R. |
| | Oswald
Registered: Apr 2002
Posts: 5023 |
if you have enough memory you can also have everything precalculated :) but there is no better way to divide on c64 thats sure.. :)
|
| | CreaMD
Registered: Dec 2001
Posts: 3035 |
Cool stuff.
Roman |
| | Krill
Registered: Apr 2002
Posts: 2851 |
the lacking accuracy of this approach is not to be ignored with 3d code. i might add that the best accuracy with this algorithm can be achieved by using the binary logarithm, but i guess you've already taken that one.
for more accuracy, the best thing is a plain binary division routine that calculates one bit of the result per loop pass, being as exact as you want it to be.
if you never coded such a thing, wait for some future issue of one of the disk mags where there are further chapters of my tutorial on maths in assembler published. i don't want to get into further detail here but i can send you sample routines or the chapter of that tutorial in german language, as i did not translate it to english yet. |
| | Scout
Registered: Dec 2002
Posts: 1568 |
Quote: the lacking accuracy of this approach is not to be ignored with 3d code. i might add that the best accuracy with this algorithm can be achieved by using the binary logarithm, but i guess you've already taken that one.
for more accuracy, the best thing is a plain binary division routine that calculates one bit of the result per loop pass, being as exact as you want it to be.
if you never coded such a thing, wait for some future issue of one of the disk mags where there are further chapters of my tutorial on maths in assembler published. i don't want to get into further detail here but i can send you sample routines or the chapter of that tutorial in german language, as i did not translate it to english yet.
Hi
I appreciate if you could send those routines/info.
It's no problem if the text is in German coz I can read it a bit (I'm dutch).
Thanx.
R.
|
| | Cruzer
Registered: Dec 2001
Posts: 1048 |
I have always precalculated everything => no need for division and other nasty realtime stuff :)
But if I'm going to make something in realtime one day, Perffs algorithm seems pretty smart... But if it's inaccurate how about this:
x/y = x*(1/y)
Then you just need a table of 1/x and to do a multiply. The latter can be done with this algorithm...
Let f(x) = x^2 / 4.
Then a*b = f(a+b) - f(a-b)
So with a table of squares/4 and 1/x you should be able to divide using adc and sbc. |
| | Nightlord
Registered: Jan 2003
Posts: 131 |
cruzer:
i have thought of division using your approach but i ran into a problem that 1/x table starts to get inacurate for large x values. besides the f(a+b) and f(a-b) also becomes troublesome. since now typical a values are decimals and typical b values are fixed point reals between 0 and 1
i have the feeling that this should be fixable by playing with the tables but couldn't figure out how? so how do you handle these issues.
thanx in advance |
| | Cruzer
Registered: Dec 2001
Posts: 1048 |
Nightlord: Do you need to divide by large numbers when doing 3d projections?
Any division routine on c64 is bound to be inaccurate or slow, so I guess you'll have to live with some inaccuracy. Then the trick is to make it just accurate enough to suit your particular need.
But as I said, I always precalculate everything, so I haven't got any experience with implementing a division routine myself.
And btw, I think Perff's routine looks smarter than mine. |
| | SIDWAVE
Account closed
Registered: Apr 2002
Posts: 2238 |
If you only want to make a 3D starfield end nothing else, just do like the game DELTA.
set a sprite and give it a different X position every rasterline.
The sprite is just made of a single vertical line in 1 color.
Do different sinus tables for each sprite, and it will look quite real to you.
|
| | Nightlord
Registered: Jan 2003
Posts: 131 |
cruzer:
at the moment i am using division for 2d line slope calculation. it is precalculated in my engine which makes it bound between 0 to $40 and i store slopes as $xx.xx so i spend $40 x $40 x 2 = $2000 bytes which hurts my mind. since i don't use eors to fill and i fill line by line (this lets me to make surfaces with dithered patterns unlike eor)i need to take one step in y everyline and left/right slope steps in left and right boundaries of line. for this reason i could not find a way to use stephen judd's add-dx-until-bigger-than-dy approach. (it isnot applicable for me when dx > dy)
can any body give some fresh ideas.
|
| | Oswald
Registered: Apr 2002
Posts: 5023 |
nightlord: my guess is, that using logs and exps is accurate enough for your routine, either u use this one, or the bloody slow, but hell accurate bit shifter killer one. (you can find this on steve judd's page)
and speaking of patterns, you can have 13 different ones if you use eorfilling a bit more trickyer: have 2 eor buffers, one with the odd, one with the even lines of your screen. then you draw each line into both of them, but with a different colour pattern for the pixelmasks. then you just eor to the screen in 2 go the odd and even lines, and voila you have patterns :) |
| | Cruzer
Registered: Dec 2001
Posts: 1048 |
When dx>dy I guess you need to do an integer division for the Bresenham algorithm, although I have never tried to implement it on the c64. Maybe you can make some LUT-stuff for this, that doesn't fill too much. |
| | White Flame
Registered: Sep 2002
Posts: 136 |
You technically don't need to do a divide.
All's that Bresenham does is keep track of the slope's fraction in a numerator/denominator pair instead of in fixed point. When the numerator gets bigger than the denominator, you bump by 1 pixel and subtract the denominator from the numerator. If you change that comparison to a "while(num>=denom) { ++x; num-=denom; }" loop (with --x for the other direction obviously), the same general code will work no matter the slope. Whether that's faster or slower than doing a divide before-hand (to figure out the integer portion of the delta so you don't need to loop the num>demon compare) depends on the slope of the line. Where the delta is large, a divide is faster. However, if the delta is large, the line won't be very long, so it sort of evens out. The code is a lot shorter and easier to manage if you do a compare loop. |
| | Nightlord
Registered: Jan 2003
Posts: 131 |
oswald: you are the master of odd/even lines
thx everyone |
| | St0fF
Registered: Oct 2002
Posts: 40 |
actually we found out some really important things here in ilmenau/germany:
the basic v2-routines are the fastest approaches for accurate calculations ... just gotta optimize those routines (or make them fixpoint)
the approach heavily beats the so highly cheered up intel math library on PCs ... just watch some stockholme syndrome releases ;) I translated the basic code to math for abductee, then he implemented it in SSE2 - fucking fast, fastest routines in the world! |
| | ready.
Registered: Feb 2003
Posts: 441 |
Don't know if anybody has come up with it, but here is my idea. I needed to do object-size-reduction in 2D-graphics, let's say a circle and let's say I want to reduce its radius to make it look fall into the screen farther and farther, smaller and smaller: I have to divide the radius by an increasing number, of course. We all know it's easy to divide by 2, 4, 8, 16, and so on, just using ROR or better LSR for shorter execution time. So it's easy to load a number into the accumulator and get 1/2, 1/4, 1/8... so I get the 50%, 25%, 12.5%....if I add these values toghether I can get 12.5%, 25%, 12.5+25=37.5%, 50%, 50+12.5=62.5%, 50+25=75%, .... of the original value.
I thought 12.5%, 25%, 50% could be enough for me, but more accurancy can be obtained doing more LSRs to get 1/16=6.25%, 1/32 and so on, and then add them toghether to get the percentage you want. Example:
lda SOMETHING
ldy #$00
sty v1
sty v2
lsr
h1: tay ;y is the 50% value of SOMETHING
lsr
h2: sta v2 ;v2 is the 25% value of SOMETHING
lsr
h3: sta v1 ;v1 is the 12.5% value of SOMETHING
lda (v1) ,y
sta v1
ldy v2
lda (v1),y ; I get the y+v2+v1=(50+25+12.5)% of SOMETHING
You have the result into the accumulator. A table should be created with numbers from 0 to $ff or even less if you need. v1 points to it.
Then you need a part which puts into the h1, h2, h3 locations the NOP code to exclude the percentage you don't want to add.
|
| | _V_
Account closed
Registered: Jan 2002
Posts: 124 |
Funny thing, division... staying out of log-space, you get Cruzer's multiplication of x*(1/y). Now, 1/y is also called the "inverse" of y, so basically, x/y = x*(1/y) is the same as saying "find the inverse of y, then scale the result by a factor of x".
An interesting fact about the inverse function [ let's call this function f: f(t) = 1/t ] is that its graph is symmetrical with respect to the primary diagonal [with point (1,1) on the diagonal]. This means that if (for some value of t) we find a point u so that
f(t) = u
then also f(u) = t
For example, f(2) = 1/2 = 0.5
but also f(0.5) = 1/0.5 = 2
[ also f(f(t)) = t as (obviously) 1/(1/t) = t ]
Anyway, this symmetry could be useful to reduce the size of the "1/y table" as described by Cruzer (or boost its accuracy by being able to add more intermediate points to the table). For example, if a 1/y table is stored as follows:
location 1/y
1 1
2 0.5
3 0.333...
4 0.25
5 0.2
etcetera,
then it's easy to find the inverse of 3: just look for the value 3 at the "location" column, see which "1/y" value is linked to it, and you get 0.333....
To the find the inverse of 0.333..., the process is reversed: look for the value 0.333... at the "1/y" column, see which "location" value is linked to it, and voila, you get 3.
In pseudo-code, you'd get something like this:
Given y, search inverse of y:
(1) Is y >= 1?
YES: goto (2)
NO : goto (3)
(2) search for y in "location" column, inverse of y is linked value in "1/y" column
(3) search for y in "1/y" column, inverse of y is linked value in "location" column
This might help Nightlord with the accuracy problem he encountered while trying Cruzer's solution. Personally, I prefer logarithms (of some base number, n) as demonstrated by Perff, as they're so elegant:
logn(x/y) = logn(x) - logn(y)
[ by the way, in Europe log(x) usually refers to base 10 (Brigg logarithm) while ln(x) refers to base e (Napier logarithm). Which got me confused about Perff's formula, but then I remembered that the US writes ln(x) as log(x) (they never heard about Brigg) =P. In what follows, I'll be using the European notation. ]
One disadvantage of logn-functions is that they don't have any easy symmetries as the one for 1/t. For large t values, ln(t) diverges very slowly to infinity, whereas it diverges to minus infinity extremely rapidly as t gets close to 0. Also logn(1) is always 0. So, you'll have to use a rather large log-table with high accuracy as t gets smaller than 1 and reduced accuracy as t gets larger and larger (the difference between ln(11000) and ln(10000) is slightly less than 0.1, for example).
There is an interesting connection between ln(t) and 1/t, by the way, which involves derivatives:
D[ln(t)] = 1/t (t>0)
The derivative of ln(t) is exactly 1/t (with t>0). So by calculating 1/t, you know exactly how fast ln(t) changes near t. Of course, derivatives aren't exactly A-material for the good old c64, so I'll leave it at this :). |
| | Graham
Account closed
Registered: Dec 2002
Posts: 990 |
Hmm nobody posted the normal binary division by now, so I'll do it:
ASL $FD
ROL
LDX #$08
.loop1
CMP $FC
BCC *+4
SBC $FC
ROL $FD
ROL
DEX
BNE .loop1
LDX #$08
.loop2
CMP $FC
BCC *+4
SBC $FC
ROL $FE
ASL
DEX
BNE .loop2
Divides the value in $FD by the value in $FC, 8 bit integer result in $FD, the first 8 fraction bits are in $FE.
Ofcourse both loops should be unrolled :) I didn't want to write down the unrolled code here. |
| | tecM0
Registered: Jan 2002
Posts: 40 |
Quote: One of the methods I have used the most when dividing (and multiplying) in assembler is using logarithm.
I use the fact that:
x/y = e^(log x - log y)
You have one table containing values of log(x), and another table of e^x
When you want to divide the content of the X-reg with the Y-reg you do:
sec
lda logtable,x
sbc logtable,y
tax
lda exptable,x
That's it! :)
(You must of course take care of overflows etc.)
One downside using this is that the result is not 100% accurate, but it's fast and you can also use it for multiplying. (do an add instead of sbc)
logtable-values stay small if i made it right...hmmm
but how you calculate the expttable?
i wanna divide 255 by 255 but precalc
in basic fails due ?exp(x) goes BIIIG fast.
T.
|
| | WVL
Registered: Mar 2002
Posts: 886 |
Quote: logtable-values stay small if i made it right...hmmm
but how you calculate the expttable?
i wanna divide 255 by 255 but precalc
in basic fails due ?exp(x) goes BIIIG fast.
T.
short answer :
don't use log10 or ln, but choose a smart base for your logarithm instead. Hint : make sure the maximum value of your log-table is $ff (I always use 2 bytes, so my max value is $ffff). This way you get the most accuracy from your table. |
| | tecM0
Registered: Jan 2002
Posts: 40 |
log (basic v2) stay under 6 for $ff isnt this right?
the problem is the EXP :(
T. |
| | WVL
Registered: Mar 2002
Posts: 886 |
Quote: log (basic v2) stay under 6 for $ff isnt this right?
the problem is the EXP :(
T.
But you want log($FF) to be $FF for biggest accuracy. Values ranging from 0->6 give no room to undo the logarithm.
To find out what your base should be, you can use the following equation :
loga(x)/loga(b)=logb(x)
as an example, we want logb(255)=255, what should b be?
logb(255)=255=log10(255)/log10(b) =>
log10(b)=log10(255)/255=0,009437 =>
b=1,02 (roughly)
So log1,02(255)=255.
in other words, for c64 stuff, forget about using 10 or e as base, 10^(log10(255)/255) is the real base of choice :D
|
| | AmiDog
Registered: Mar 2003
Posts: 97 |
Has anyone tried the "multiply with scaled reciprocal" approach for division?
Using an 8-bit table, some bitshifting and a 8x8 multiply (using the table of squares method) one should get a 100% accurate quotient. |
| | John West
Account closed
Registered: Aug 2004
Posts: 4 |
Quoting AmiDogHas anyone tried the "multiply with scaled reciprocal" approach for division?
It's not quite 100% accurate, but it's never more than one out. It gets more results exactly right if the reciprocal table is rounded, so you're calculating floor((x*floor(0.5+256/y))/256). The rounding can be built into the table, so it's no more work. |
| | Graham
Account closed
Registered: Dec 2002
Posts: 990 |
1 off is ok, remember that you have an error in that range anyway due to rounding.
Oh and btw, I use X*(1/Y) too when I need a fast div. |
| | John West
Account closed
Registered: Aug 2004
Posts: 4 |
... and if you don't mind it being up to two out, you don't need to calculate the low byte of the multiply. That saves a few cycles. |
| | AmiDog
Registered: Mar 2003
Posts: 97 |
Well, you do need an 9-bit reciprocal to get a 100% accurate result, but as the msb of the reciprocal is always one, you don't need to store it anywhere, and a multiply by 256 is just an add. Some C-code:
/* 9-bit reciprocal table, but msb is always one and not stored */
unsigned char r0[128];
void build_table(void)
{
unsigned long i, d, r;
for (i=0; i<128; i++) {
d = i | 0x80; /* 1 < d < 2 */
r = 0x80000000 / (d << 8); /* 1 / d -> 1/2 < r < 1 */
r += 0x80; /* round the reciprocal */
r <<= 1; /* remove constant 1-bit */
r0[i] = (unsigned char)(r >> 8);
}
}
unsigned char div8(unsigned char a, unsigned char b)
{
int log = 7;
/* division by one need to be handled separately */
if (b == 1) {
return a;
}
/* normalize b so that 1 <= b < 2 */
while (!(b & 0x80)) {
b <<= 1;
log--;
}
/* if b is a power of two */
if (b == 0x80) {
return a >> log;
}
/* msb is always one, so don't use it for indexing (save some memory) */
unsigned char r = r0[b & 0x7f];
/* now multiply */
unsigned short ar = (unsigned short)a * (unsigned short)r;
ar >>= 8;
ar += (unsigned short)a; /* ar is now an 9-bit value */
ar >>= (log + 1);
return (unsigned char)ar;
}
This has been tested and is 100% accurate. I don't know how fast/slow it will be when translated to 6502 code. Maybe adding a table or two for the normalizing step could help. |
| | AmiDog
Registered: Mar 2003
Posts: 97 |
Just noticed that a divide by zero will cause an endless loop, so better replace:
if (b == 1) {
return a;
}
with
if (b <= 1) {
return a;
}
to avoid the endless loop. The result of a divide by zero is usually considered undefined anyway, so it doesn't really matter what we return... |
| | AmiDog
Registered: Mar 2003
Posts: 97 |
I just couldn't control myself, here's some 6502 code:
_divu_8
lda div_b
cmp #2
bcs + ; >= 2
lda div_a
rts
+ ldx #8
- dex
asl
bcc -
bne +
lda div_a
- lsr
dex
bne -
rts
+ tay
lda r0_table,y
ldy div_a
sta zp8_1
sta zp8_2
eor #$ff
sta zp8_3
sta zp8_4
sec
lda (zp8_1),y
sbc (zp8_3),y
lda (zp8_2),y
sbc (zp8_4),y
clc
adc div_a
ror
- lsr
dex
bne -
rts
div_a
.byte $0
div_b
.byte $0
r0_table
.byte $01,$00,$fd,$00,$f9,$00,$f5,$00,$f1,$00,$ed,$00,$ea,$00,$e6,$00
.byte $e2,$00,$df,$00,$db,$00,$d8,$00,$d5,$00,$d1,$00,$ce,$00,$cb,$00
.byte $c8,$00,$c4,$00,$c1,$00,$be,$00,$bb,$00,$b8,$00,$b5,$00,$b3,$00
.byte $b0,$00,$ad,$00,$aa,$00,$a7,$00,$a5,$00,$a2,$00,$9f,$00,$9d,$00
.byte $9a,$00,$98,$00,$95,$00,$93,$00,$90,$00,$8e,$00,$8b,$00,$89,$00
.byte $87,$00,$84,$00,$82,$00,$80,$00,$7e,$00,$7b,$00,$79,$00,$77,$00
.byte $75,$00,$73,$00,$71,$00,$6f,$00,$6d,$00,$6b,$00,$69,$00,$67,$00
.byte $65,$00,$63,$00,$61,$00,$5f,$00,$5d,$00,$5b,$00,$59,$00,$58,$00
.byte $56,$00,$54,$00,$52,$00,$51,$00,$4f,$00,$4d,$00,$4b,$00,$4a,$00
.byte $48,$00,$47,$00,$45,$00,$43,$00,$42,$00,$40,$00,$3f,$00,$3d,$00
.byte $3c,$00,$3a,$00,$39,$00,$37,$00,$36,$00,$34,$00,$33,$00,$31,$00
.byte $30,$00,$2f,$00,$2d,$00,$2c,$00,$2a,$00,$29,$00,$28,$00,$26,$00
.byte $25,$00,$24,$00,$22,$00,$21,$00,$20,$00,$1f,$00,$1d,$00,$1c,$00
.byte $1b,$00,$1a,$00,$19,$00,$17,$00,$16,$00,$15,$00,$14,$00,$13,$00
.byte $12,$00,$10,$00,$0f,$00,$0e,$00,$0d,$00,$0c,$00,$0b,$00,$0a,$00
.byte $09,$00,$08,$00,$07,$00,$06,$00,$05,$00,$04,$00,$03,$00,$02,$00
This will divide two 8-bit numbers in some 90-150 cycles. The code can easily be extended to handle larger dividends. |
| | AmiDog
Registered: Mar 2003
Posts: 97 |
Let me bore you with an optimized version:
; divide acc by y, result in acc
_divu_8
ldx t0_table,y
stx b1+1
ldx t1_table,y
beq +
ldy r0_table,x
sta zp8_1
sta zp8_2
eor #$ff
sta zp8_3
sta zp8_4
sec
lda (zp8_1),y
sbc (zp8_3),y
lda (zp8_2),y
sbc (zp8_4),y
clc
adc zp8_1
ror
+ sec
b1 bcs b1
lsr
lsr
lsr
lsr
lsr
lsr
lsr
rts
.align $100
r0_table
.byte $01,$00,$fd,$00,$f9,$00,$f5,$00,$f1,$00,$ed,$00,$ea,$00,$e6,$00
.byte $e2,$00,$df,$00,$db,$00,$d8,$00,$d5,$00,$d1,$00,$ce,$00,$cb,$00
.byte $c8,$00,$c4,$00,$c1,$00,$be,$00,$bb,$00,$b8,$00,$b5,$00,$b3,$00
.byte $b0,$00,$ad,$00,$aa,$00,$a7,$00,$a5,$00,$a2,$00,$9f,$00,$9d,$00
.byte $9a,$00,$98,$00,$95,$00,$93,$00,$90,$00,$8e,$00,$8b,$00,$89,$00
.byte $87,$00,$84,$00,$82,$00,$80,$00,$7e,$00,$7b,$00,$79,$00,$77,$00
.byte $75,$00,$73,$00,$71,$00,$6f,$00,$6d,$00,$6b,$00,$69,$00,$67,$00
.byte $65,$00,$63,$00,$61,$00,$5f,$00,$5d,$00,$5b,$00,$59,$00,$58,$00
.byte $56,$00,$54,$00,$52,$00,$51,$00,$4f,$00,$4d,$00,$4b,$00,$4a,$00
.byte $48,$00,$47,$00,$45,$00,$43,$00,$42,$00,$40,$00,$3f,$00,$3d,$00
.byte $3c,$00,$3a,$00,$39,$00,$37,$00,$36,$00,$34,$00,$33,$00,$31,$00
.byte $30,$00,$2f,$00,$2d,$00,$2c,$00,$2a,$00,$29,$00,$28,$00,$26,$00
.byte $25,$00,$24,$00,$22,$00,$21,$00,$20,$00,$1f,$00,$1d,$00,$1c,$00
.byte $1b,$00,$1a,$00,$19,$00,$17,$00,$16,$00,$15,$00,$14,$00,$13,$00
.byte $12,$00,$10,$00,$0f,$00,$0e,$00,$0d,$00,$0c,$00,$0b,$00,$0a,$00
.byte $09,$00,$08,$00,$07,$00,$06,$00,$05,$00,$04,$00,$03,$00,$02,$00
t0_table
.fill $100,0
t1_table
.fill $100,0
_divu_8_setup
ldy #1
next
tya
ldx #$ff
- inx
asl
bcc -
sta t1_table,y
txa
sta t0_table,y
iny
bne next
rts
This will divide two 8-bit numbers in 42-106 cycles (including the jsr used to call this).
I'll stop now, I promise ;) |
| | Krill
Registered: Apr 2002
Posts: 2851 |
Yes, and now please re-read some former posts by other people for somewhat better approaches :) |
| | AmiDog
Registered: Mar 2003
Posts: 97 |
Well, it depends on how you define "better". If you need accuracy, this method will give you that (all bits are correct) in 40-100 cycles (I've made a few more optimizations). Compare this with the unrolled binary division in a post above which will use 244-274 cycles...
Also, you could easily extend this code to perform a 16bit/8bit division, or even 32bit/8bit division by just extending the multiplication and the bitshifting at the end, the reciprocal and shiftcount stays the same. Somehow I doubt you'll be able to do that using log-tables :) |
| | Style
Registered: Jun 2004
Posts: 498 |
pwned krill
:) |
| | Krill
Registered: Apr 2002
Posts: 2851 |
Ok, i'm sorry, after a closer investigation, the code is much more elegant and better than i suspected after a quick peek. :) Well done. |
|
