[CSDb] - User Forums - Fast large multiplies

Welcome to our latest new user Jedfox ! (Registered 2024-05-28)

You are not logged in - nap

CSDb User Forums

Forums > C64 Coding > Fast large multiplies

2012-06-09 19:45

Repose

Registered: Oct 2010
Posts: 222

Fast large multiplies

I've discovered some interesting optimizations for multiplying large numbers, if the multiply routine time depends on the bits of the mulitplier. Usually if there's a 1 bit in the multiplier, with a standard shift and add routine, there's a "bit" more time or that bit.
The method uses several ways of transforming the input to have less 1 bits. Normally, if every value appears equally, you average half 1 bits. In my case, that becomes the worst case, and there's about a quarter 1 bits. This can speed up any routine, even the one that happens to be in rom, by using pre- and post- processing of results. The improvement is about 20%.
Another speedup is optimizing the same multiplier applied to multiple multiplicands. This saves a little in processing the multiplier bits once. This can save another 15%.
Using the square table method will be faster but use a lot of data and a lot of code.
Would anyone be interested in this?

... 144 posts hidden. Click here to view all posts....

2017-04-08 04:43

Repose

Registered: Oct 2010
Posts: 222

Post correction would be slower on the 16x16, not a long enough run of add/sub.

2017-04-08 05:00

Repose

Registered: Oct 2010
Posts: 222

Partials cheatsheet

                        y1  y0
                        x1  x0
                        ------
                 x0*y0h x0*y0l
          x1*y0h x1*y0l
          x0*y1h x0*y1l
   x1*y1h x1*y1l

----------------
24x24bits

                    x2     x1    x0
                    y2     y1    y0
                    ---------------
                        y0x0h y0x0l
                  y0x1h y0x1l
            y0x2h y0x2l
                  y1x0h y1x0l
            y1x1h y1x1l
      y1x2h y1x2l
            y2x0h y2x0l
      y2x1h y2x1l
y2x2h y2x2l

These facts are useful to estimating the time of any size calc:

Number of columns is 2*n, n is bytes of each number.
Rows of additions is like 1 3 3 1, 1 3 5 5 3 1, 1 3 5 7 7 5 3 1 for 16,24 and 32 bit mults, and each one being f(x)-g(x), so really double that number of add-subs.
The total add-subs is n^2*2. (each is about 10 cycles).
Number of times to change (or pointers to set) of the multiplier is n, then each one is used >n times with the multiplicand (ldy multiplicand), when doing in column order (tbd total).

2017-04-08 05:44

Repose

Registered: Oct 2010
Posts: 222

Changes in multiplicand is 2*n (in my case, ldy y(x) ).
eg ldy y0
...
ldy y1
...
ldy y0
...
etc.

2017-04-08 07:42

ChristopherJam

Registered: Aug 2004
Posts: 1381

I've replaced all the subtractions with additions by using g(x)=$4000-(x*x/4) and offsetting my start number, and rolled the carry correction in to the addition sequence.

Removing the CLCs is probably not worthwhile for 32x32, as there are only 33 of them, so I'd have to spend considerably less than 4 cycles per output byte on the post correction (unlikely).

I'm down to around 800 cycles for a 32x32 now, 776 cycles for zero times zero.

2017-04-08 07:43

ChristopherJam

Registered: Aug 2004
Posts: 1381

fo=open("tables.inc","w")
lo=lambda x:x&255
hi=lambda x:(x>>8)
f=lambda x:x*x//4
g=lambda x:(0x4000-f(x-255))&0xffff
dumpArrayToA65(fo, "flo", [lo(f(i)) for i in range(512)])
dumpArrayToA65(fo, "fhi", [hi(f(i)) for i in range(512)])
dumpArrayToA65(fo, "glo", [lo(g(i)) for i in range(512)])
dumpArrayToA65(fo, "ghi", [hi(g(i)) for i in range(512)])
dumpArrayToA65(fo, "id",  [lo(  i ) for i in range(512)])

fo=open("mc.inc","w")
mAcc=0
for i in range(4):
    for j in range(4):
        mAcc-=0x40<<(8*(1+i+j))

initialValue = [((mAcc>>s)&0xff) for s in range(0,64,8)]

def addB(yv,zp,tb):
    global lasty
    if yv!=lasty:
        print("""    ldy mT2+{yv}""".format(yv=yv), file=fo)
        lasty=yv
    print("""    adc ({zp}),y""".format(zp=zp), file=fo)
    if tb<7:
        print("""    bcc *+4:inx:clc""", file=fo)
    else:
        print("""    bcc *+3:clc""", file=fo)

lasty=None
for tb in range(8):
    print(""" ; tb={tb} """.format(tb=tb),file=fo)
    if tb==0:
        print("""    clc """,file=fo)
        print("""    ldx#0 """,file=fo)
        print("""    lda #${iv:02x} """.format(iv=initialValue[tb]),file=fo)
    else:
        print("""    txa""", file=fo)
        if tb<7:
            print("""    ldx#0 """,file=fo)
        print("""    adc#${iv:02x}""".format(iv=initialValue[tb]), file=fo)

        if initialValue[tb]>0xef:
            print("""    bcc *+4:inx:clc""", file=fo)


    for j in range(4):
        i=tb-j
        if i in [0,1,2,3]:
            addB(i, "zp_fl{j}".format(j=j), tb)
            addB(i, "zp_gl{j}".format(j=j), tb)
        i=tb-j-1
        if i in [0,1,2,3]:
            addB(i, "zp_fh{j}".format(j=j), tb)
            addB(i, "zp_gh{j}".format(j=j), tb)
    print("""    sta mRes+{tb}""".format(tb=tb), file=fo)

fo.close()

2017-04-08 07:48

ChristopherJam

Registered: Aug 2004
Posts: 1381

(obviously also need four sets of

    lda mT1+ 0
    sta zp_fl0
    sta zp_fh0
    eor#255
    sta zp_gl0
    sta zp_gh0

every time the multiplier changes (included in my cycle times above),
plus also some init code to set the high bytes of the zero page pointers before a set of multiplications is performed (not included in my timings))

2017-04-08 08:37

Repose

Registered: Oct 2010
Posts: 222

Good job, that's right in the range of what I thought was possible.

I have an improvement; instead of trashing A to change the multiplier, you can prestuff pointers with the 4 multipliers.

by offset $4000, doesn't that reduce the domain?

Correction is fast, it's only

stx z3
sec
sbc z3
sta z2

Also yours shouldn't be any faster than my approach from what I can tell, though I do have some ideas to speed up adds again.. we'll see :)

2017-04-08 09:50

ChristopherJam

Registered: Aug 2004
Posts: 1381

Quoting Repose

Good job, that's right in the range of what I thought was possible.

Thanks!

Quote:

I have an improvement; instead of trashing A to change the multiplier, you can prestuff pointers with the 4 multipliers.

Doing :)

Quote:

by offset $4000, doesn't that reduce the domain?

Nah, second table only contains x**2/4 for x in -255 to 255, so it already maxed out at $3f80

Quote:

Correction is fast, it's only

stx z3 sec sbc z3 sta z2

True, but that's an extra 64 cycles, and removing the CLCs saves at most 64 cycles, sometimes as little as zero (if the branches skip over them all)

Quote:

Also yours shouldn't be any faster than my approach from what I can tell, though I do have some ideas to speed up adds again.. we'll see :)

Yes, there should be an equivalent that mixes ADC and SBC, I just found it easier to wrap my brain around the edge cases and carry handling by converting it to ADC only. I'll be interested to see what you come up with.

2017-04-08 10:19

ChristopherJam

Registered: Aug 2004
Posts: 1381

Oh, faster correction:

  sec
  sbc id,x
  sta z2

Still not gonna do it, mind ;)

2017-04-08 10:59

Repose

Registered: Oct 2010
Posts: 222

About the correction, I think you're adding things up wrong. I only use correction for those columns where it's faster, and I found the break even at 7 adds, so it should work. All but the outer 1 or 2 columns can use it.

Let's say in the middle columns where there's 14 adds per column, that's 28 cycles half the time saved from not using CLC, or 14 cycles on average, vs 8 cycles for correction, it still saves 6.

I actually found the stats for the carries, most of them are about half, but adding a higher proportion of high bytes gives less carries.