home *** CD-ROM | disk | FTP | other *** search
/ Fujiology Archive / fujiology_archive_v1_0.iso / !FALCON / LINEOUT / OUT.ZIP / SOURCE.ZIP / MATH.I < prev    next >
Text File  |  2002-11-03  |  17KB  |  740 lines
  1. ; Math libs. all stuff in here, sin/cos/tan/exp/atan/random,etc.
  2.  
  3. ;======= equates ===========================================================
  4.  
  5. * >WARNING< for these equ's: when using a new 'sintbllen' you must
  6. * recalculate 'cos1' and 'sin1'!
  7. sintbllen:    equ    2048        * MUST BE A EXPONENTIAL VALUE OF 2!
  8. sin1:        equ    13176774    * sin(2π/2048)*2^32
  9. cos1:        equ    4294947083    * cos(2π/2048)*2^32
  10.  
  11. Math.EXP_ITERATIONS:    =    10
  12.  
  13. ; This is really awful.. To avoid distortion near angles of +-pi/4, you
  14. ; need at least 20 iterations!
  15. Math.SATAN:        =    1            ; fast but 'Sneaky ATAN'
  16. Math.ATAN_ITERATIONS:    =    40
  17. Math.RAD2DEG:        =    sintbllen*5000/31415    ;degs/2pi
  18.  
  19. Calc_NextRandom:    MACRO
  20.     move.l    random,d0
  21.     move.l    d0,d1
  22.     mulu.w    d0,d0
  23.     eor.l    d1,d0
  24.     addq.l    #7,d0
  25.     move.l    d0,random
  26.     ENDM
  27.  
  28. getNextRandom:
  29.     Calc_NextRandom
  30.     rts
  31.  
  32. * Very fast and accurate squareroot algorithm.
  33. * Quite lengthy, though: 66 bytes.
  34. * INPUT: d1.l: value to calculate the squareroot of (integer)
  35. * OUTPUT: d0.l: squareroot of value (16.16 fixed point)
  36. Math.sqrt:
  37. CALC_ATARISQRT:
  38.     moveq    #1,d2
  39.     ror.l    #2,d2
  40.     moveq    #$F,d3
  41. .loop1:    cmp.l    d2,d1
  42.     bgt.s    .endloop1
  43.     add.l    d1,d1
  44.     lsr.l    #1,d2
  45.     dbf    d3,.loop1
  46.     moveq    #0,d0
  47.     bra.s    .is_null
  48. .endloop1:
  49.  
  50.     sub.l    d2,d1
  51.     move.l    d2,d0
  52.     lsr.l    #1,d2
  53. .loop2:    lsr.l    #1,d2
  54.     add.l    d2,d0
  55.     cmp.l    d0,d1
  56.     bgt.s    .endloop2
  57.     sub.l    d2,d0
  58.     add.l    d1,d1
  59.     dbf    d3,.loop2
  60.     bra.s    .end
  61. .endloop2:
  62.  
  63.     sub.l    d0,d1
  64.     add.l    d2,d0
  65.     add.l    d1,d1
  66.     dbf    d3,.loop2
  67.  
  68. .end:    add.l    d0,d0
  69.     addi.l    #$00008000,d0
  70. .is_null:
  71.     rts
  72.  
  73. *==========================================================================
  74. * Sinewave table generator.
  75. * By EarX/~fUn~, 10-5-1998
  76. * 68020 or higher is required!
  77. *==========================================================================
  78.  
  79. * Macro that returns the modulo of a given angle.
  80. * INPUT: angle: type: data-register (word) or RAM (word)
  81. Do_SinModulo:    MACRO    angle
  82.         andi.w    #sintbllen-1,\1
  83.         ENDM
  84.  
  85. * Macro that returns sine & cosine of a given angle.
  86. * PRECONDITION: INIT_SINETABLE has been called!
  87. * INPUT: base: type: address-register or address or relative address
  88. *        inpreg: type: data-register or address-register (lower word)
  89. *                contains: angle (0=0°, sintbllen=360°)
  90. * OUTPUT: sinreg: type: data-register (long) or address-register
  91. *                 contains: sine value (signed: -32768 to 32767)
  92. *         cosreg: type: data-register (long) or address-register
  93. *                 contains: cosine value (signed: -32768 to 32767)
  94. Get_SinCos:    MACRO    base,inpreg,sinreg,cosreg
  95.         movem.w    (\1,\2.w*4),\3/\4
  96.         ENDM
  97.  
  98. * Macro that returns sine of a given angle.
  99. * PRECONDITION: INIT_SINETABLE has been called!
  100. * INPUT: base: type: address-register or address or relative address
  101. *        inpreg: type: data-register or address-register (lower word)
  102. *                contains: angle (0=0°, sintbllen=360°)
  103. * OUTPUT: sinreg: type: data-register (long) or address-register
  104. *                 contains: sine value (signed: -32768 to 32767)
  105. Get_Sin:    MACRO    base,inpreg,sinreg
  106.         move.w    (\1,\2.w*4),\3
  107.         ENDM
  108.  
  109. * Macro that returns cosine of a given angle.
  110. * PRECONDITION: INIT_SINETABLE has been called!
  111. * INPUT: base: type: address-register or address or relative address
  112. *        inpreg: type: data-register or address-register (lower word)
  113. *                contains: angle (0=0°, sintbllen=360°)
  114. * OUTPUT: cosreg: type: data-register (long) or address-register
  115. *                 contains: cosine value (signed: -32768 to 32767)
  116. Get_Cos:    MACRO    base,inpreg,cosreg
  117.         move.w    2(\1,\2.w*4),\3
  118.         ENDM
  119.  
  120. * Creates the a combined sine and cosine table for quick fetching.
  121. * Macro is exactly 96 bytes in length :-)
  122. * INPUT: a0: address of sine_tbl
  123. Init_SineTable:    MACRO
  124.         moveq    #$ffffffff,d0        * /cos(0)=1
  125.         lsr.l    #1,d0            * \(=$7fffffff)
  126.         moveq    #0,d1            * sin(0)=0
  127.         move.l    #sin1,d6
  128.         move.w    #sintbllen/4-1,d7
  129.  
  130. .genlop:    swap    d0            * Get high-word of cosa
  131.         swap    d1            * Get high-word of sina
  132.         move.w    d1,2+(sintbllen)*3(a0)    * Copy sina in cos-4th quadrant
  133.         move.w    d0,sintbllen*1(a0)    * Copy cosa in sin-2nd quadrant
  134.         sub.w    d1,2+(sintbllen)*1(a0)    * Copy -sina in cos-2nd quadrant
  135.         sub.w    d0,sintbllen*3(a0)    * Copy -cosa in sin-4th quadrant
  136.         sub.w    d0,2+(sintbllen)*2(a0)    * Copy -cosa in cos-3rd quadrant
  137.         sub.w    d1,sintbllen*2(a0)    * Copy -sina in sin-3rd quadrant
  138.         move.w    d1,(a0)+        * Save sina (16 bit signed value) in first quadrant
  139.         move.w    d0,(a0)+        * Save cosa (16 bit signed value) in first quadrant
  140.         swap    d0            * Change cosa back to fixedpoint
  141.         swap    d1            * Change sina back to fixedpoint
  142.         move.l    d1,d4            * / Backup sina 
  143.         move.l    d0,d5            * | and cosa
  144.         move.l    d1,d2            * | for use in
  145.         move.l    d0,d3            * \ multiplications.
  146.         mulu.l    d6,d3:d1        * d3:=sin1*sina
  147.         mulu.l    #cos1,d2:d0        * d2:=cos1*cosa
  148.         mulu.l    d6,d1:d5        * d0:=sin1*cosa
  149.         mulu.l    #cos1,d0:d4        * d1:=cos1*sina
  150.         sub.l    d3,d2            * d2:=(cos1*cosa)-(sin1*sina)
  151.         add.l    d0,d1            * sina:=(sin1*cosa)+(cos1*sina)
  152.         move.l    d2,d0            * cosa:=(cos1*cosa)-(sin1*sina)
  153.         dbra    d7,.genlop
  154.         ENDM
  155.  
  156. * Creates a tangens table by using the sine/cosine table.
  157. * INPUT: sincos_tbl: address register
  158. * OUTPUT: tan_tbl: address register
  159. Init_TanTable:    MACRO    sincos_tbl,tan_tbl
  160.         move.w    #sintbllen-1,d7
  161. .gentanloop:    movem.w    (/1)+,d0-d1
  162.         move.w    d0,d0
  163.         beq.s    .notandiv
  164.         lsl.l    #8,d1
  165.         divs.w    d0,d1
  166. .notandiv:    move.w    d1,(/2)+
  167.         dbra    d7,.gentanloop
  168.         ENDM
  169.  
  170. ;------------------------------------------------------------------------------
  171. ; exp(x) using Taylor series:
  172. ;
  173. ; exp(x)=1+x+(1/2!)x^2+(1/3!)x^3+..
  174. ; exp(x)=Sum[ n : 0<=n : x^n/n! ]
  175. ;
  176. ;------------------------------------------------------------------------------
  177. ; Proof of correctness by variation of constant N:
  178. ;
  179. ; exp(x,N)= Sum[ n : 0<=n<=N : x^n/n! ]
  180. ;
  181. ; exp(x,N+1)= exp(x,N) + x^(N+1)/(N+1)!
  182. ;
  183. ;                x      x^N
  184. ; = exp(x,N) + ----- * -----
  185. ;               N+1     N!
  186. ;
  187. ; => { T.n = x^n/n! }
  188. ;
  189. ;------------------------------------------------------------------------------
  190. ; Main invariants:
  191. ;
  192. ; exp(x,0)=T.0=1
  193. ; T.(N+1)=x/(N+1)*T.n
  194. ; exp(x,N+1) = exp(x,N) + x/(N+1)*T.n
  195. ;
  196. ;------------------------------------------------------------------------------
  197.  
  198. ; Unsigned 16b fractional implementation. Inputs and outputs fixedpoint.
  199. ; Reasonably accurate, at least: when |x|>1 => |x|<#iterations*2!!
  200. ; AND when numbers can be extremely large (even more than doubleprec!).
  201. ; |x|<=1 always works fast and accurate.
  202. ; Needs 020 or higher..
  203. ;
  204. ; Better idea is to precalc e and 1/e and put them in tables.
  205. ; Then only power function is needed. This is simple cos of the short
  206. ; fractional taylor approx and the highly optimisable multiplier with
  207. ; precalced tables!! Maybe next time..
  208. ;
  209. ; INPUT:
  210. ; d0.l=x (16:16)
  211. ; OUTPUT:
  212. ; d1.l=exp(x) (16:16)
  213. Math.exp:
  214. ; todo: power function and table stuff!
  215.  
  216. ; Good fractional exp function. Even works with few iterations.
  217. Math.fracExp:
  218.     IFNE    1
  219.  
  220.     bsr    Fix32.toFloat
  221. ; d0.l=x (f)
  222.     moveq    #1,d7
  223.     move.l    #$7FFFFF00,d6            ; d6.l=T.0=1.0 (f)
  224.     move.l    d0,d5                ; d5.l=x (f)
  225.     movea.l    d6,a0                ; a0=exp(x,0) (f)
  226.  
  227. .loop:    move.l    d6,d1                ; d1.l=T.n (f)
  228.     move.l    d5,d0                ; d0.l=x (f)
  229. ; d0.l=x (f), d1.l=T.n (f)
  230.     bsr    Float.mul
  231.     move.l    d0,d6
  232. ; d6.l=x*T.n (f)
  233.     move.l    d7,d0
  234.     bsr    Int.toFloat
  235.     move.l    d0,d1
  236.     move.l    d6,d0
  237. ; d0.l=x*T.n (f), d1.l=n+1 (f)
  238.     bsr    Float.div    
  239. ; d0.l=T.(n+1)=x/(n+1)*T.n (f)
  240.     move.l    d0,d6                ; Store T.(n+1)
  241.     move.l    a0,d1
  242. ; d1.l=exp(x,n)
  243.     bsr    Float.add
  244.     movea.l    d0,a0
  245. ; d5.l=exp(x,n+1)
  246.     addq.w    #1,d7                ; n:=n+1
  247.     cmpi.w    #Math.EXP_ITERATIONS,d7
  248.     blt.s    .loop
  249.  
  250.     bsr    Float.toFix32
  251. ; d0.l=exp(x) (16:16)
  252.  
  253.     ELSE
  254.  
  255.     moveq    #1,d7
  256.     move.l    #$00010000,d1
  257.     move.l    d1,d2
  258.  
  259. .loop:    muls.l    d0,d3:d2
  260.     move.w    d3,d2
  261.     swap    d2
  262. ; d2.l=x*T.n (16:16)
  263.     divs.l    d7,d2
  264. ; d2.l=T.(n+1)=x/(N+1)*T.n
  265.     add.l    d2,d1
  266. ; d1.l=exp(x,n+1)=exp(x,n)+T.(n+1)
  267.     addq.w    #1,d7
  268.     cmpi.w    #Math.EXP_ITERATIONS,d7
  269.     blt.s    .loop
  270.  
  271.     ENDC
  272.     rts
  273.  
  274. ;------------------------------------------------------------------------------
  275. ; Definition of arctan(x) using Taylor series:
  276. ;
  277. ; arctan(x)=x-x^3/3+x^5/5-..., (-1<x<+1)
  278. ;
  279. ; arctan(x)= Sum[ n : 0<=n : (-1)^n*x^(1+2n)/(1+2n) ], (-1<x<+1)
  280. ;
  281. ;------------------------------------------------------------------------------
  282. ; Proof of correctness by variation of constant N:
  283. ;
  284. ; arctan(x,N)= Sum[ n : 0<=n<=N : (-1)^n*x^(1+2n)/(1+2n) ]
  285. ;
  286. ; arctan(x,N+1)= arctan(x,N) + (-1)^(N+1)*x^(1+2(N+1))/(1+2(N+1))
  287. ;
  288. ; p.n=1+2n
  289. ; p.0=1
  290. ; p.(n+1)=p.n+2
  291. ;
  292. ; T.n=x^p.n
  293. ; T.0=x^p.0=x^1=x
  294. ; T.(n+1)=x^p.(n+1)=x^(p.n+2)=x^p.n*x*x=T.n*x*x
  295. ;
  296. ;------------------------------------------------------------------------------
  297. ; Main invariants:
  298. ;
  299. ; arctan(x,0)= (-1)^0*x^1/1= T.0= x
  300. ;
  301. ; arctan(x,N+1)= arctan(x,N) + (-1)^(N+1)*T.N*x*x/(p.N+2)
  302. ;------------------------------------------------------------------------------
  303.  
  304. ; Unsigned 16b fractional implementation. Inputs and outputs fixedpoint.
  305. ; Reasonably accurate and works on old 68000.
  306. ; INPUT:
  307. ; d0.l=x (16:16) (-256.0<x<+256.0)
  308. ; OUTPUT:
  309. ; d1.l=arctan(x) (radians) (16:16)
  310. Frac.arctan:
  311.  
  312.     move.l    d0,-(sp)
  313.     bpl.s    .sign_okay
  314.     neg.l    d0
  315. .sign_okay:
  316.     cmpi.l    #$00010000,d0
  317.     bne.s    .calc_quadrant
  318.  
  319. ; Output pi/4 (arctan(1)=pi/4).
  320.     move.l    #31415<<14/10000,d1
  321.     tst.l    (sp)+
  322.     bpl.s    .sign_okay2
  323.     neg.l    d1
  324. .sign_okay2:
  325.     rts
  326.  
  327. .calc_quadrant:
  328.     bgt.s    .adjust
  329.     bsr.s    .calc
  330.     tst.l    (sp)+
  331.     bpl.s    .end_sign
  332.     neg.l    d1
  333. .end_sign:
  334.     rts
  335.  
  336. .adjust:
  337.     move.l    #$01000000,d1
  338.     lsr.l    #8,d0
  339.     divu.w    d0,d1
  340.     move.w    d1,d0                ; d0.w=1/x (frac)
  341. ;    moveq    #$FFFFFFFF,d1
  342. ;    divu.l    d0,d1
  343. ;    move.w    d1,d0
  344.  
  345.     bsr.s    .calc
  346.     move.l    #31415<<15/10000,d2
  347.     sub.l    d1,d2
  348.     move.l    d2,d1                ; d1.l=pi/2-arctan(1/x)
  349.     tst.l    (sp)+
  350.     bpl.s    .end_sign2
  351.     neg.l    d1
  352. .end_sign2:
  353.     rts
  354.     
  355. ; d0.l=x (16:16, 0<=x<1)
  356. ; OUTPUT:
  357. ; d1.l=arctan(x) (16:16, non-negative)
  358. .calc:
  359.  
  360.     IFNE    Math.SATAN
  361.  
  362. ; atan(a) =~ x/(1+0.28*x*x)
  363.     move.w    d0,d2
  364.     mulu.w    d2,d2
  365.     swap    d2                ; d2.w=x*x (frac)
  366.     mulu.w    #65536*28/100,d2        ; d2.l=0.28*x*x (frac)
  367.     move.w    #1,d2
  368.     swap    d2                ; d2.l=1+0.28*x*x (16:16)
  369.     swap    d0                ; d0.l=x<<16 (16:16)
  370.     divu.l    d2,d0                ; d0.l=x/(1+0.28*x*x)
  371.     move.l    d0,d1
  372.     
  373.     ELSE
  374.  
  375.     clr.l    d1
  376.     move.w    d0,d1                ; d1.l=arctan(x,n)={n=0}=x (frac)
  377.     moveq    #1,d2                ; d2.l=p.n=2n+1
  378.     move.w    d0,d3                ; d3.l=T.n={n=0}=x (frac)
  379.     mulu.w    d0,d0
  380.     swap    d0                ; d0.w=x*x (frac)
  381.     moveq    #Math.ATAN_ITERATIONS/2-1,d6    ; d6.w=counter
  382.  
  383. .loop:
  384. ; Calculate p.(n+1).
  385.     addq.w    #2,d2                ; d2.l=p.n+2=p.(n+1)
  386.  
  387. ; Calculate T.(n+1).
  388.     mulu.w    d0,d3                ; d3.l=T.(n+1)=x*x*T.n (16:16)
  389.     swap    d3                ; d3.w=T.(n+1)=x*x*T.n (frac)
  390.     clr.l    d5
  391.     move.w    d3,d5
  392.     divu.w    d2,d5                ; d5.w=x*x*T.n/(p.n+2) (frac)
  393.  
  394. ; Calculate arctan(x,n+1).
  395.     sub.w    d5,d1
  396.  
  397. ; Calculate p.(n+1).
  398.     addq.w    #2,d2                ; d2.l=p.n+2=p.(n+1)
  399.  
  400. ; Calculate T.(n+1).
  401.     mulu.w    d0,d3                ; d3.l=T.(n+1)=x*x*T.n (16:16)
  402.     swap    d3                ; d3.w=T.(n+1)=x*x*T.n (frac)
  403.     clr.l    d5
  404.     move.w    d3,d5
  405.     divu.w    d2,d5                ; d5.w=x*x*T.n/(p.n+2) (frac)
  406.  
  407. ; Calculate arctan(x,n+1).
  408.     add.w    d5,d1
  409.     dbra    d6,.loop
  410.  
  411.     ENDC
  412.  
  413.     rts
  414.  
  415. ; Uses two input system..
  416. ; INPUT:
  417. ; d0.l=x (int)
  418. ; d1.l=y (int)
  419. Frac.atan2:
  420.     clr.l    d3
  421.     cmp.l    d0,d1
  422.     bhi.s    .sorted
  423.     bne.s    .go_on
  424.     moveq    #$FFFFFFFF,d0
  425.     bra.s    .a_calced
  426. .go_on:    exg    d0,d1
  427.     not.l    d3
  428. .sorted:swap    d0
  429.     divu.w    d1,d0                ; d0.w=a=x/y (frac)
  430. .a_calced:
  431.  
  432. ; atan(a) =~ a/(1+0.28*a*a)
  433.     move.w    d0,d2
  434.     mulu.w    d2,d2
  435.     swap    d2                ; d2.w=a*a (frac)
  436.     mulu.w    #65536*28/100,d2        ; d2.l=0.28*a*a (frac)
  437.     move.w    #1,d2
  438.     swap    d2                ; d2.l=1+0.28*a*a (16:16)
  439.     swap    d0                ; d0.l=a<<16 (16:16)
  440.     divu.l    d2,d0                ; d0.l=a/(1+0.28*a*a)
  441.     move.l    d0,d1
  442.  
  443.     tst.l    d3
  444.     beq.s    .end
  445.     move.l    #31415<<15/10000,d2
  446.     sub.l    d1,d2
  447.     move.l    d2,d1                ; d1.l=pi/2-arctan(1/a)
  448.  
  449. .end:    rts
  450.  
  451.  
  452.     IFNE    0
  453.  
  454. ; Floating point implementation, accuracy overkill, you all know I'm sick.
  455. ; Needs 68020 or higher.
  456. ;
  457. ; integers: n, p.n
  458. ; reals:    x, x*x, T.n, arctan(x,n)
  459. ;
  460. ; INPUT:
  461. ; d0.l=x (float)
  462. ; OUTPUT:
  463. ; d1.l=arctan(x) (radians) (float)
  464. Float.arctan:
  465.     cmpi.l    #$40000001,d0
  466.     bne.s    .calc_quadrant
  467.  
  468. ; Output pi/4 (arctan(1)=pi/4).
  469.     move.l    #$6487ED00,d0
  470.     rts
  471.  
  472. .calc_quadrant:
  473.     move.l    #$40000001,d1
  474.     move.l    d0,-(sp)
  475.     bsr    Float.sub            ; d0.l=x-1.0 (float)
  476.     move.l    d0,d1                ; d1.l=x-1.0 (float)
  477.     move.l    (sp)+,d0            ; d0.l=x (float)
  478.     tst.l    d1
  479.     bgt.s    .adjust
  480.     bra.s    .calc
  481.  
  482. .adjust:move.l    d0,d1                ; d1.l=x (float)
  483.     move.l    #$40000001,d0            ; d0.l=1.0 (float)
  484.     bsr    Float.div            ; d0.l=1/x (float)
  485.     bsr.s    .calc                ; d1.l=arctan(1/x) (float)
  486.     move.l    #$6487ED01,d0            ; d0.l pi/2 (float)
  487.     bsr    Float.sub
  488.     move.l    d0,d1                ; d1.l=pi/2-arctan(1/x) (float)
  489.     rts
  490.     
  491. .calc:    move.l    d0,d1
  492.     move.l    d0,-(sp)
  493.     bsr    Float.mul
  494.     move.l    d0,d4                ; d4.l=x*x (float)
  495.     move.l    (sp)+,d0
  496.     sub.l    a6,a6                ; a6=n=0
  497.     moveq    #1,d7                ; d7.l=p.n=2n+1
  498.     move.l    d0,d1                ; d1.l=arctan(x,n)={n=0}=T.n (float)
  499.     move.l    d0,d3                ; d3.l=T.n (float)
  500.  
  501. .loop:    addq.w    #2,d7                ; d7.l=(p.n+2)=p.(n+1) (int)
  502.  
  503. ; Calculate T.(n+1).
  504.     move.l    d1,-(sp)
  505.     move.l    d3,d0
  506.     move.l    d4,d1
  507.     bsr    Float.mul
  508. ; d0.l=T.(n+1)=x*x*T.n (float)
  509.     move.l    d0,-(sp)
  510.     move.l    d7,d0
  511.     bsr    Int.toFloat
  512.     move.l    d0,d1
  513.     move.l    (sp),d0
  514. ; d0.l=T.(n+1)=x*x*T.n (float)
  515. ; d1.l=p.(n+1) (float)
  516.     bsr    Float.div
  517.     move.l    (sp)+,d3
  518. ; d3.l=T.(n+1)=x*x*T.n (float)
  519. ; d0.l=T.(n+1)=x*x*T.n/p.(n+1) (float)
  520.     move.l    (sp)+,d1
  521.  
  522. ; Calculate (-1)^(n+1)*T.(n+1).
  523. .sign:    addq    #1,a6
  524.     move.l    a6,d5
  525.     btst    #0,d5
  526.     beq.s    .end_sign
  527.     move.l    d1,-(sp)
  528.     bsr    Float.neg
  529.     move.l    (sp)+,d1
  530. .end_sign:
  531.  
  532. ; Calculate arctan(x,n+1).
  533. ; d0.l=(-1)^(n+1)*T.(n+1) (float)
  534. ; d1.l=arctan(x,n) (float)
  535.     movem.l    d3-d4,-(sp)
  536.     bsr    Float.add
  537.     movem.l    (sp)+,d3-d4
  538.     move.l    d0,d1
  539. ; d1.l=arctan(x,n+1) (float)
  540.  
  541.     cmpa.w    #Math.ATAN_ITERATIONS,a6
  542.     blt.s    .loop
  543.     rts
  544.  
  545.     ENDC
  546.  
  547. ; INPUT:
  548. ; d0.l=src1 (float)
  549. ; d1.l=src2 (float)
  550. ; OUTPUT:
  551. ; d0.l=dest (float) src1*src2
  552. Float.mul:
  553.     clr.l    d2
  554.     clr.l    d3
  555.     move.b    d0,d2                ; d2.b=exp src1
  556.     move.b    d1,d3                ; d3.b=exp src2
  557.     sub.l    d2,d0                ; d0.l=mantisse src1
  558.     sub.l    d3,d1                ; d1.l=mantisse src2
  559.     muls.l    d0,d0:d1
  560.     add.l    d1,d1
  561.     addx.l    d0,d0                ; d0.l=mantisse of src1*src2
  562.     subq.b    #1,d2
  563.     add.l    d1,d1
  564.     addx.l    d0,d0
  565.     bvc.s    .end_normalization
  566.     addx.b    d0,d0
  567.     ror.l    d0
  568.     addq.b    #1,d2
  569. .end_normalization:
  570.     add.b    d3,d2                ; d2.b=exp of src1*src2
  571.     move.b    d2,d0
  572.     rts
  573.  
  574. ; INPUT:
  575. ; d0.l=divident (float)
  576. ; d1.l=divisor (float)
  577. ; OUTPUT:
  578. ; d0.l=quotient (float) divident/divisor
  579. Float.div:
  580.     clr.l    d2
  581.     clr.l    d3
  582.     move.b    d0,d2                ; d2.b=exp src1
  583.     move.b    d1,d3                ; d3.b=exp src2
  584.     sub.l    d2,d0                ; d0.l=mantisse src1
  585.     sub.l    d3,d1                ; d1.l=mantisse src2
  586.     sub.b    d3,d2                ; d2.b=exp of src1*src2
  587.     asr.l    #2,d0                ; divident<divisor, cos of result limit.
  588.     clr.l    d3
  589.     divs.l    d1,d0:d3            ; d0.l=mantisse of src1*src2
  590.     move.l    d3,d0
  591.     add.l    d0,d0
  592.     bvc.s    .end_normalization
  593.     move.l    d3,d0
  594.     addq.b    #1,d2
  595. .end_normalization:
  596.     move.b    d2,d0
  597.     rts
  598.  
  599. ; INPUT:
  600. ; d0.l=src1 (24b mantisse:8b exponent)
  601. ; d1.l=src2 (24b mantisse:8b exponent)
  602. ; OUTPUT:
  603. ; d0.l=dest (24b mantisse:8b exponent) src1-src2
  604. Float.sub:
  605.     move.b    d1,d3
  606.     clr.b    d1
  607.     neg.l    d1
  608.     move.b    d3,d1
  609.  
  610. ;
  611. ; m: mantissa, p: power
  612. ;
  613. ; a+b~= { m=m(a)+m(b)>>[p(a)-p(b)], p=p(a); p(a)>=p(b)
  614. ;       { m=m(b)+m(a)>>[p(b)-p(a)], p=p(b); p(a)<p(b)
  615. ;
  616. ; INPUT:
  617. ; d0.l=src1 (24b mantisse:8b exponent)
  618. ; d1.l=src2 (24b mantisse:8b exponent)
  619. ; OUTPUT:
  620. ; d0.l=dest (24b mantisse:8b exponent) src1+src2
  621. Float.add:
  622.     clr.l    d2
  623.     clr.l    d3
  624.     move.b    d0,d2                ; d2.b=p(a)
  625.     move.b    d1,d3                ; d3.b=p(b)
  626.     sub.l    d2,d0                ; d0.l=m(a)
  627.     sub.l    d3,d1                ; d1.l=m(b)
  628.     cmp.b    d3,d2                ; p(a)>=p(b)??
  629.     bge.s    .end_exg
  630.     exg.l    d0,d1
  631.     exg.l    d2,d3
  632. .end_exg:
  633.     move.b    d2,d4
  634.     sub.b    d3,d2                ; d2.b=p(a)-p(b)
  635.     asr.l    d2,d1                ; d1.l=m(b)>>[p(a)-p(b)]
  636.     addq.b    #1,d4
  637.     add.l    d1,d0                ; d0.l=m(a)+m(b)>>[p(a)-p(b)]
  638.     bvc.s    .normalize_loop
  639.     addx.b    d0,d0                ; if carried..
  640.     ror.l    #1,d0                ; ..make negative
  641.     move.b    d4,d0
  642.     rts
  643.  
  644. .normalize_loop:
  645.     subq.b    #1,d4
  646.     move.l    d0,d2
  647.     add.l    d0,d0
  648.     bvc.s    .normalize_loop
  649. .end_normalization:
  650.     move.l    d2,d0
  651.     move.b    d4,d0                ; d0.l={m(a)+m(b)>>[p(a)-p(b)],p=p(a)}=a+b
  652.     rts
  653.  
  654. ; INPUT:
  655. ; d0.l=src (24b mantisse:8b exponent)
  656. ; OUTPUT:
  657. ; d0.l=-src (24b mantisse:8b exponent)
  658. Float.neg:
  659.     move.b    d0,d1
  660.     clr.b    d0
  661.     neg.l    d0
  662.     move.b    d1,d0
  663.     rts
  664.  
  665. ; INPUT:
  666. ; d0.l=src (24m:8e)
  667. ; OUTPUT:
  668. ; d0.l=src (32b signed)
  669. Float.toInt:
  670.     move.b    d0,d1
  671.     clr.b    d0
  672.     subi.b    #31,d1                ; d1.b=exp-[log(r)+7]
  673.     bpl.s    .limit_0
  674.     neg.b    d1
  675.     asr.l    d1,d0
  676.     rts
  677. .limit_0:
  678.     clr.l    d0
  679.     rts
  680.  
  681. ; INPUT:
  682. ; d0.l=src (24m:8e)
  683. ; OUTPUT:
  684. ; d0.l=src (16:16 sgn)
  685. Float.toFix32:
  686.     move.b    d0,d1
  687.     clr.b    d0
  688.     subi.b    #15,d1                ; d1.b=exp-[log(r)+7]+16
  689.     bpl.s    .limit_0
  690.     neg.b    d1
  691.     asr.l    d1,d0
  692.     rts
  693. .limit_0:
  694.     clr.l    d0
  695.     rts
  696.  
  697. ; INPUT:
  698. ; d0.l=src (32b signed)
  699. ; OUTPUT:
  700. ; d0.l=src (24m:8e)
  701. Int.toFloat:
  702.     moveq    #31,d1
  703.     tst.l    d0
  704.     beq.s    .limit_0
  705. .normalize_loop:
  706.     move.l    d0,d2
  707.     add.l    d0,d0
  708.     dbvs    d1,.normalize_loop
  709.     move.l    d2,d0
  710.     move.b    d1,d0                ; d0.l={m(a)+m(b)>>[p(a)-p(b)],p=p(a)}=a+b
  711.     rts
  712. .limit_0:
  713.     move.l    #$40000080,d0    
  714.     rts
  715.  
  716. ; Converts 16:16 fixed point number to floating point.
  717. ; INPUT:
  718. ; d0.l=src (32b signed)
  719. ; OUTPUT:
  720. ; d0.l=src (24m:8e)
  721. Fix32.toFloat:
  722.     moveq    #31,d1
  723.     tst.l    d0
  724.     beq.s    .limit_0
  725. .normalize_loop:
  726.     move.l    d0,d2
  727.     add.l    d0,d0
  728.     dbvs    d1,.normalize_loop
  729.     move.l    d2,d0
  730.     move.b    d1,d0                ; d0.l={m(a)+m(b)>>[p(a)-p(b)],p=p(a)}=a+b
  731.     subi.b    #16,d0                ; /65536
  732.     rts
  733. .limit_0:
  734.     move.l    #$40000080,d0    
  735.     rts
  736.  
  737.     BSS
  738.  
  739. random:    DS.L    1
  740.