home *** CD-ROM | disk | FTP | other *** search
/ Kosovo Orphans' Appeal Charity CD / KosovoOrphansAppeal.iso / archimedesworld_cd2 / acornanswers / _scape2 / s / subrouts1 < prev   
Text File  |  1996-02-21  |  19KB  |  754 lines
  1.  
  2. max16    *    &7fffffff    ;max positive number in 16 bit fixed point format (2^31-1)/65536
  3. min16    *    &80000000    ;min negative number                  -2^31   /65536
  4. one    *    65536
  5.  
  6.  
  7.  
  8. ; rbbcinc
  9. ; a leaf APCS function
  10. ;
  11. ; C prototype:
  12. ; int rbbcinc(int r, int k)
  13. ;
  14. ; given limits 0 <= r < 2^k, return (r ñ 1) where ñ denotes a reversed bit ordering increment,
  15. ; subject to the stated limits.
  16. ; NB in this case, for the function f: n -> {for (k=c=0; c<n; c++, k=rbbcinc(k, l); return k;},
  17. ; we have f(f(n))=n.
  18. ; NB2 algo requires 2 <= k <= 32, but doesn't check for this - BEWARE!
  19. ; (for k=1 get no action taken, while for any other bad k get code executed at unintended address,
  20. ;  hence unpredictable & likely system fatal).
  21. ;
  22.  
  23.         EXPORT  rbbcinc
  24.  
  25. rbnsta  DCB     "rbbcinc", 0
  26.         ALIGN
  27. rbnend  DCD     &ff000000 + rbnend - rbnsta
  28.  
  29. rbbcinc
  30.  
  31.         rbbc    a1, a2, a3
  32.         MOVS    pc, lr
  33.  
  34.  
  35.  
  36. ; pow16
  37. ; a leaf APCS function
  38. ;
  39. ; C prototype:
  40. ; int pow16(int a, int b)
  41. ;
  42. ; returns 65536 * (a/65536)^(b/65536)
  43. ;
  44. ; is about 5 to 40 times faster than FPEmulator working with double floats
  45. ;
  46. ; nb if a<0 require b an integer
  47. ;
  48. ; nb2 unexpected return values:
  49. ; a=b=0                returns 1    actual value ill defined
  50. ; a<0, b non integer        returns 0    actual value complex
  51. ; abs(bln(a)) > max16        unknown errors likely (to avoid, restrict b to ▒(2^11)one, roughly)
  52. ; a=0, b<0            returns max16    actual value +infinity
  53. ; a & b st result out of range    unknown result
  54. ;
  55. ; nb precise nature of errors for other values not yet confirmed acceptable, though should be okay except
  56. ; perhaps in extreme cases (as of 22/11/95)
  57. ;
  58. ; Few subsequent notes on error:
  59. ;
  60. ; For result much bigger than one:
  61. ; since exp16 only yields 16 significant bits, pow16 gives no more than this, hence for result much bigger
  62. ; than one get increasing error, eg for a around 64one & b=1.03125one get errors st answer (around 73one)
  63. ; is only correct to top 15 or 16 bits (ie error upto around one/1024 to one/256)
  64. ;
  65. ; For 0<=a<=one, & b eg 2one, 3one, 13.125one etc get error typically no more than in low 1 or 2 bits
  66. ;
  67. ; For a large (say a>one) & b<0 get error typically in only lowest bit or no error
  68. ;
  69. ;   ****  Strongly recommend this fn is only used where accuracy not crucial or with limited range ****
  70. ;   ****  of arguments where you can confirm yourself accuracy in that range is adequate.          ****
  71. ;
  72.  
  73.         EXPORT  pow16
  74.  
  75. pwnsta  DCB     "pow16", 0
  76.         ALIGN
  77. pwnend  DCD     &ff000000 + pwnend - pwnsta
  78.  
  79. pow16
  80.  
  81.     CMP    a2, #0
  82.     MOVEQ    a1, #&10000        ;if b=0 return one directly
  83.     MOVEQS    pc, lr
  84.     CMP    a2, #&10000
  85.     MOVEQS    pc, lr            ;handle trival case b=one directly, by returning a
  86.     CMP    a1, #0
  87.     MOVGT    ip, #0
  88.     BGT    pow16_apos
  89.     BNE    pow16_aneg
  90.     CMP    a2, #0
  91.     MOV    a1, #0            ;if a=0 and b>0 return 0
  92.     MOVLT    a1, #max16+1        ;if a=0 and b<0 return max16
  93.     SUBLT    a1, a1, #1
  94.     MOVS    pc, lr
  95. pow16_aneg
  96.     MOVS    ip, a2, LSL #16
  97.     MOVNE    a1, #0            ;if a<0 and b not an integer return 0
  98.     MOVNES    pc, lr
  99.     AND    ip, a2, #&10000        ;if a<0 then go on to evaluate (65536 * (-a/65536)^(b/65536))
  100.     RSB    a1, a1, #0        ;with sign subsequently forced to + or - according to whether b/one
  101. pow16_apos                ;is even or odd ( eg (-2.5)^3 is just (-1)^3 * 2.5^3 )
  102.     STMFD    sp!, {v1, v2, lr}
  103.     MOV    v1, ip            ;now evaluate (65536 * (a/65536)^(b/65536)), for a>0
  104.     MOV    v2, a2            ;via exp16( ln16(a) * b / 65536 )
  105.     BL    ln16            ;nb this uses property that exp(b.log(a)) = exp(log(a^b)) = a^b
  106.         mul16    a1, v2, a1, a2, a3, a4
  107.     BL    exp16
  108.     CMP    v1, #0
  109.     RSBNE    a1, a1, #0        ;finally switch sign on answer if originally a<0 and b/one was odd
  110.     LDMFD    sp!, {v1, v2, pc}^
  111.  
  112.  
  113.  
  114. ; ln16
  115. ; a leaf APCS function
  116. ;
  117. ; C prototype:
  118. ; int ln16(int a)
  119. ;
  120. ; returns 65536 * ln (a/65536)
  121. ;
  122. ; is about 30 times faster than FPEmulator working with double floats
  123. ;
  124.  
  125.         EXPORT  ln16
  126.  
  127. lnnsta  DCB     "ln16", 0
  128.         ALIGN
  129. lnnend  DCD     &ff000000 + lnnend - lnnsta
  130.  
  131. ln16
  132.  
  133.     CMP    a1, #0            ;if a<=0 return min16
  134.     MOVLE    a1, #min16
  135.     MOVLES    pc, lr            ;else most significant bit is between b30 & b0
  136.     GBLA    counter
  137. counter    SETA    30            ;find msb & store (msb_index - 15) in a4
  138.     WHILE    counter > 1
  139.     TST    a1, #1<<counter
  140.     MOVNE    a4, #counter-15
  141.     BNE    ln16_gotmsb
  142. counter    SETA    counter-1
  143.     WEND
  144.     TST    a1, #1<<1
  145.     MOVNE    a4, #1-15
  146.     MOVEQ    a4, #0-15
  147. ln16_gotmsb                ;if we set z = a1/(2^a4), will have ln(a/one) = ln(z/one * 2^a4)
  148.     CMP    a4, #2            ;                          = ln(z/one) + a4*ln2
  149.     SUBGT    a2, a4, #2        ;with z in range [0.5,1)
  150.     MOVGT    a1, a1, LSR a2        ;so calc a1 = 4*( a1/(2^a4) ) - 3*one
  151.     RSBLE    a2, a4, #2        ;which puts a1 in range ▒one suitable for approximation by poly
  152.     MOVLE    a1, a1, LSL a2        ;leaving us to calculate a4*one*ln2 + one*ln((a1+3one)/4one)
  153.     SUB    a1, a1, #&30000
  154.     ADD    a2, a1, a1, LSL #2    ;start polynomial approximation calculation on a1
  155.     RSB    a2, a2, a1, LSL #10
  156.     MOV    a2, a2, ASR #16        ;for details of how this works see comments against similar code
  157.     SUB    a2, a2, #&000e00    ;in exp16 function, directly below this function
  158.     SUB    a2, a2, #&000034
  159.         MOV    ip, a1
  160.     mul16c    a2, ip, a2, a3
  161.     ADD    a2, a2, #&003200
  162.     ADD    a2, a2, #&000031
  163.         MOV    ip, a1
  164.     mul16c    a2, ip, a2, a3
  165.     SUB    a2, a2, #&00e200
  166.     SUB    a2, a2, #&0000f2
  167.         MOV    ip, a1
  168.     mul16c    a2, ip, a2, a3
  169.     ADD    a2, a2, #&050000
  170.     ADD    a2, a2, #&005500
  171.     ADD    a2, a2, #&000065
  172.         MOV    ip, a1
  173.     mul16c    a2, ip, a2, a3
  174.     SUB    a2, a2, #&040000
  175.     SUB    a2, a2, #&009a00    ;end of polynomial approximation calculation
  176.     SUB    a1, a2, #&000059    ;a1 now holds (2^20)*(approximated value)
  177.     ADD    ip, a4, a4, LSL #1
  178.     ADD    a2, ip, a4, LSL #3    ;now add in the a4*one*ln2 term using an explicit binary expansion
  179.     ADD    a2, a4, a2, LSL #4    ;of approximately one*ln2
  180.     RSB    a3, a4, a4, LSL #3
  181.     ADD    a2, a3, a2, LSL #4    ;nb we actually add approx (2^20)*ln2, & then divide by 16
  182.     ADD    a2, a4, a2, LSL #3    ;to reduce truncation error
  183.     ADD    a1, a1, a2, LSL #5
  184.     ADD    a1, a1, ip, ASR #1
  185.     MOV    a1, a1, ASR #4
  186.         MOVS    pc, lr
  187.  
  188.  
  189.  
  190. ; exp16
  191. ; a leaf APCS function
  192. ;
  193. ; C prototype:
  194. ; int exp16(int a)
  195. ;
  196. ; returns 65536 * exp (a/65536)
  197. ;
  198. ; is about 45 times faster than FPEmulator working with double floats
  199. ;
  200.  
  201.         EXPORT  exp16
  202.  
  203. epnsta  DCB     "exp16", 0
  204.         ALIGN
  205. epnend  DCD     &ff000000 + epnend - epnsta
  206.  
  207. exp16                    ;this fn returns a value with only about 17 significant binary
  208.                     ;digits - eg for 'a' small or negative, get near maximal accuracy
  209.                     ;but for 'a' large (upto around 10one) where exp has value
  210.     CMP    a1, #12*65536        ;~20000one can get error upto roughly 0.1one
  211.     MOVGT    a1, #max16+1
  212.     SUBGT    a1, a1, #1
  213.     MOVGTS    pc, lr
  214.     CMP    a1, #-12*65536
  215.     MOVLT    a1, #0
  216.     MOVLTS    pc, lr            ;otherwise know |a| <= 12one
  217.     RSB    a2, a1, a1, LSL #3    ;now set a = a/ln2 using explicit mul by 2_1.01110001010101000111
  218.     ADD    a3, a1, a1, LSL #2    ;nb only need 1st 20 digits given range on a
  219.     ADD    a3, a3, a3, LSL #4    ;nb2 this calc is approximate only, though will usually yield an
  220.     ADD    a1, a2, a1, LSL #4    ;answer correct to the stored 16 binary places (else error one/65536)
  221.     ADD    a1, a1, a3, ASR #10
  222.     ADD    a1, a1, a2, ASR #16    ;thus exp(original a)  = exp(new a * ln2) = 2^(new a), eval'd below:
  223.     ADD    a1, a1, #8
  224.     MOV    a1, a1, ASR #4        ;calc a = a/ln2 complete
  225.     MOV    a4, a1, ASR #16        ;a4 now holds integer component of a
  226.     BIC    a1, a1, a4, LSL #16    ;a1 now holds fractional component in range [0,1)*one
  227.     CMP    a4, #15            ;do another check on a to ensure exp(a) is in range
  228.     MOVGE    a1, #max16+1        ;& if not return maximal or minimal value as appropriate
  229.     SUBGE    a1, a1, #1
  230.     MOVGE    pc, lr
  231.     CMP    a4, #-16
  232.     MOVLT    a1, #0
  233.     MOVLTS    pc, lr            ;else need to return value (one * 2^(a4 + a1/one))
  234.     MOV    a1, a1, LSL #1        ;             = (one * 2^(a1/one) * 2^a4)
  235.     SUB    a1, a1, #&10000        ;convert a1 to lie in [-1,-1), then apply poly approximation:
  236.     RSB    a2, a1, a1, LSL #3
  237.     ADD    a2, a1, a2, LSL #6
  238.     MOV    a2, a2, ASR #15        ;a2 = 0.0008561*(2^20)*(a1/one)
  239.     ADD    a2, a2, #&002800
  240.     ADD    a2, a2, #&00007e    ;a2 = 0.0008561*(2^20)*(a1/one) + 0.0098857*(2^20)
  241.         MOV    ip, a1
  242.     mul16c    a2, ip, a2, a3        ;a2 = 0.0008561*(2^20)*(a1/one)^2 + 0.0098857*(2^20)*(a1/one)
  243.     ADD    a2, a2, #&015000
  244.     ADD    a2, a2, #&000be0    ; etc
  245.         MOV    ip, a1
  246.     mul16c    a2, ip, a2, a3
  247.     ADD    a2, a2, #&070000
  248.     ADD    a2, a2, #&00d700
  249.     ADD    a2, a2, #&00007e
  250.         MOV    ip, a1
  251.     mul16c    a2, ip, a2, a3        ;until we have
  252.     ADD    a2, a2, #&160000    ;a2 = (2^20) ( 0.0008561(a1/one)^4 + 0.0098857(a1/one)^3 +
  253.     ADD    a2, a2, #&00a000    ;           0.0849301(a1/one)^2 + 0.4901106(a1/one) +
  254.     ADD    a2, a2, #&0000a7    ;           1.14142138(a1/one)                           )
  255.     MOV    a1, a2, ASR #4        ;now divide by 16 removing most truncation error from above calcs,
  256.     CMP    a4, #0            ;leaving a2 = (2^16)*(approximated value)
  257.     RSBLT    a4, a4, #0
  258.     MOVLT    a1, a1, ASR a4        ;finally carry out the (a2 = a2 * 2^a4) step
  259.     MOVLTS    pc, lr
  260.     MOVS    a1, a1, LSL a4
  261.     MOVMI    a1, #max16+1
  262.     SUBMI    a1, a1, #1
  263.         MOVS    pc, lr
  264.  
  265.  
  266.  
  267. ; cos16
  268. ; a leaf APCS function
  269. ;
  270. ; C prototype:
  271. ; int cos16(int a)
  272. ;
  273. ; returns 65536 * cos ((a/65536)(pi/2))
  274. ;
  275. ; is about 44 times faster than FPEmulator working with double floats
  276. ;
  277.  
  278.         EXPORT  cos16
  279.  
  280. csnsta  DCB     "cos16", 0
  281.         ALIGN
  282. csnend  DCD     &ff000000 + csnend - csnsta
  283.  
  284. cos16
  285.  
  286.     EOR    a4, a1, a1, LSL #1    ;since cos is periodic and cos((a/65536)(pi/2)) has period 4*one
  287.     TST    a4, #&20000        ;eval is simpler than for ln or exp:
  288.     MOV    a4, #0            ;we need only to truncate a to [0,4one) and then approx the function
  289.     MOVNE    a4, #1            ;via a poly
  290.     TST    a1, #&10000
  291.     MOV    a1, a1, LSL #16        ;in fact further symmetries in cos over this range allow us to
  292.     MOV    a1, a1, LSR #16        ;actually only approximate function over range [0,one)
  293.     RSBEQ    a1, a1, #&10000        ;and reconstruct it over remainder of range [one,4one) from that part
  294.     MOV    a1, a1, LSL #1        ;eg  cos( (a/65536)(pi/2) ) for a in [2one, 3one)
  295.     SUB    a1, a1, #&10000        ; = -cos( ((a-2one)/65536)(pi/2) )
  296.     RSB    a2, a1, a1, LSL #3
  297.     ADD    a2, a1, a2, LSL #8
  298.     MOV    a2, a2, ASR #16
  299.     ADD    a2, a2, #&002c00
  300.     ADD    a2, a2, #&000086
  301.         MOV    ip, a1
  302.     mul16c    a2, ip, a2, a3
  303.     SUB    a2, a2, #&00e900
  304.     SUB    a2, a2, #&0000bd
  305.         MOV    ip, a1
  306.     mul16c    a2, ip, a2, a3
  307.     SUB    a2, a2, #&030000
  308.     SUB    a2, a2, #&007c00
  309.     SUB    a2, a2, #&0000c6
  310.         MOV    ip, a1
  311.     mul16c    a2, ip, a2, a3
  312.     ADD    a2, a2, #&080000
  313.     ADD    a2, a2, #&00E200
  314.     ADD    a2, a2, #&0000BD
  315.         MOV    ip, a1
  316.     mul16c    a2, ip, a2, a3
  317.     ADD    a2, a2, #&0b0000
  318.     ADD    a2, a2, #&005000
  319.     ADD    a2, a2, #&000050
  320.     MOV    a1, a2, ASR #4
  321.     TST    a4, #1
  322.     RSBNE    a1, a1, #0
  323.         MOVS    pc, lr
  324.  
  325.  
  326.  
  327. ; sin16
  328. ; a leaf APCS function
  329. ;
  330. ; C prototype:
  331. ; int sin16(int a)
  332. ;
  333. ; returns 65536 * sin ((a/65536)(pi/2))
  334. ;
  335. ; is about 44 times faster than FPEmulator working with double floats
  336. ;
  337.  
  338.         EXPORT  sin16
  339.  
  340. snnsta  DCB     "sin16", 0
  341.         ALIGN
  342. snnend  DCD     &ff000000 + snnend - snnsta
  343.  
  344. sin16
  345.  
  346.     TST    a1, #&20000        ;eval of sin is almost identical to that of cos
  347.     MOV    a4, #0
  348.     MOVNE    a4, #1
  349.     TST    a1, #&10000
  350.     MOV    a1, a1, LSL #16
  351.     MOV    a1, a1, LSR #16
  352.     RSBNE    a1, a1, #&10000
  353.     MOV    a1, a1, LSL #1
  354.     SUB    a1, a1, #&10000
  355.     RSB    a2, a1, a1, LSL #3
  356.     ADD    a2, a1, a2, LSL #8
  357.     MOV    a2, a2, ASR #16
  358.     ADD    a2, a2, #&002c00
  359.     ADD    a2, a2, #&000086
  360.         MOV    ip, a1
  361.     mul16c    a2, ip, a2, a3
  362.     SUB    a2, a2, #&00e900
  363.     SUB    a2, a2, #&0000bd
  364.         MOV    ip, a1
  365.     mul16c    a2, ip, a2, a3
  366.     SUB    a2, a2, #&030000
  367.     SUB    a2, a2, #&007c00
  368.     SUB    a2, a2, #&0000c6
  369.         MOV    ip, a1
  370.     mul16c    a2, ip, a2, a3
  371.     ADD    a2, a2, #&080000
  372.     ADD    a2, a2, #&00E200
  373.     ADD    a2, a2, #&0000BD
  374.         MOV    ip, a1
  375.     mul16c    a2, ip, a2, a3
  376.     ADD    a2, a2, #&0b0000
  377.     ADD    a2, a2, #&005000
  378.     ADD    a2, a2, #&000050
  379.     MOV    a1, a2, ASR #4
  380.     TST    a4, #1
  381.     RSBNE    a1, a1, #0
  382.         MOVS    pc, lr
  383.  
  384.  
  385.  
  386. ; acs16
  387. ; a leaf APCS function
  388. ;
  389. ; C prototype:
  390. ; int acs16(int a)
  391. ;
  392. ; returns 65536 * 2/pi * arccos (a/65536)
  393. ;
  394. ; for a out of range (ie abs a > 65536), reset a to ▒65536 as appropriate before eval
  395. ;
  396. ; is about 24 times faster than FPEmulator working with double floats
  397. ;
  398.  
  399.         EXPORT  acs16
  400.  
  401. acnsta  DCB     "acs16", 0
  402.         ALIGN
  403. acnend  DCD     &ff000000 + acnend - acnsta
  404.  
  405. acs16
  406.  
  407.     CMP    a1, #65536
  408.     MOVGE    a1, #0
  409.     MOVGES    pc, lr
  410.     CMP    a1, #-65536
  411.     MOVLE    a1, #2*65536
  412.     MOVLES    pc, lr
  413.     STMFD    sp!, {lr}        ;now have arg in range (-one,one)
  414.     CMP    a1, #0
  415.     RSBMI    a1, a1, #0
  416.     MOVMI    a4, #1            ;b0 in a4 set iff final result r needs replacing by 2*one-r
  417.     MOVPL    a4, #0            ;now have arg in range [0,one)
  418.     SUB    a2, a1, #&b500
  419.     SUBS    a2, a2, #&0005
  420.     BLT    acs16_arglow        ;if arg < one/sqrt2 go and apply poly else first transmute into here
  421.     ORR    a4, a4, #2        ;b1 in a4 clear iff penultimate result r need be replaced by one-r
  422.     MUL    a2, a1, a1
  423.     MOV    a1, a2, LSR #16
  424.     RSB    a1, a1, #65536
  425.     sqrt16    a1, a1, a2, a3, ip, lr    ;arg transmuted (ie replaced by one-sqrt(arg^2))
  426. acs16_arglow
  427.     ADD    a3, a1, a1, ASL #2
  428.     ADD    a2, a3, a1, ASL #4
  429.     ADD    a2, a2, a3, ASL #5
  430.     ADD    a2, a2, a3, ASR    #8
  431.     ADD    a2, a2, #32
  432.     MOV    a1, a2, ASR #6        ;range transform applied (ie replaced by arg*2sqrt2-one)
  433.     SUB    a1, a1, #&10000        ;so arg lies in range [-1,1]
  434.     ADD    a2, a1, a1, LSL #2    ;now apply poly
  435.     ADD    a2, a2, a2, LSL #5
  436.     MOV    a2, a2, ASR #14
  437.     ADD    a2, a2, #&000500
  438.     ADD    a2, a2, #&0000ae
  439.         MOV    ip, a1
  440.     mul16c    a2, ip, a2, a3
  441.     ADD    a2, a2, #&000800
  442.     ADD    a2, a2, #&000088
  443.         MOV    ip, a1
  444.     mul16c    a2, ip, a2, a3
  445.     ADD    a2, a2, #&002000
  446.     ADD    a2, a2, #&00009a
  447.         MOV    ip, a1
  448.     mul16c    a2, ip, a2, a3
  449.     ADD    a2, a2, #&004600
  450.     ADD    a2, a2, #&00009a
  451.         MOV    ip, a1
  452.     mul16c    a2, ip, a2, a3
  453.     ADD    a2, a2, #&030000
  454.     ADD    a2, a2, #&00d900
  455.     ADD    a2, a2, #&0000b4
  456.         MOV    ip, a1
  457.     mul16c    a2, ip, a2, a3
  458.     ADD    a2, a2, #&030000
  459.     ADD    a2, a2, #&00ae00
  460.     ADD    a2, a2, #&000052
  461.     MOV    a1, a2, ASR #4
  462.     TST    a4, #2
  463.     RSBEQ    a1, a1, #65536
  464.     TST    a4, #1
  465.     RSBNE    a1, a1, #2*65536
  466.     LDMFD    sp!, {pc}^
  467.  
  468.  
  469.  
  470. ; asn16
  471. ; a leaf APCS function
  472. ;
  473. ; C prototype:
  474. ; int asn16(int a)
  475. ;
  476. ; returns 65536 * 2/pi * arcsin (a/65536)
  477. ;
  478. ; for a out of range (ie abs a > 65536), reset a to ▒65536 as appropriate before eval
  479. ;
  480. ; is about 24 times faster than FPEmulator working with double floats
  481. ;
  482.  
  483.         EXPORT  asn16
  484.  
  485. asnsta  DCB     "asn16", 0
  486.         ALIGN
  487. asnend  DCD     &ff000000 + asnend - asnsta
  488.  
  489. asn16
  490.  
  491.     CMP    a1, #65536
  492.     MOVGE    a1, #65536
  493.     MOVGES    pc, lr
  494.     ADDS    a2, a1, #65536        ;equiv to CMP a1,#-65536
  495.     SUBLE    a1, a1, a2        ;equiv to MOVLE a1,#-65536 (which wouldn't work due to bad constant)
  496.     MOVLES    pc, lr
  497.     STMFD    sp!, {lr}        ;now have arg in range (-one,one)
  498.     CMP    a1, #0
  499.     RSBMI    a1, a1, #0
  500.     MOVMI    a4, #1            ;b0 in a4 set iff final result needs negating
  501.     MOVPL    a4, #0            ;now have arg in range [0,one)
  502.     SUB    a2, a1, #&b500
  503.     SUBS    a2, a2, #&0005
  504.     BLT    asn16_arglow        ;if arg < one/sqrt2 go and apply poly else first transmute into here
  505.     ORR    a4, a4, #2        ;b1 in a4 set iff penultimate result r needs to be replaced by one-r
  506.     MUL    a2, a1, a1
  507.     MOV    a1, a2, LSR #16
  508.     RSB    a1, a1, #65536
  509.     sqrt16    a1, a1, a2, a3, ip, lr    ;arg transmuted (ie replaced by one-sqrt(arg^2))
  510. asn16_arglow
  511.     ADD    a3, a1, a1, ASL #2
  512.     ADD    a2, a3, a1, ASL #4
  513.     ADD    a2, a2, a3, ASL #5
  514.     ADD    a2, a2, a3, ASR    #8
  515.     ADD    a2, a2, #32
  516.     MOV    a1, a2, ASR #6        ;range transform applied (ie replaced by arg*2sqrt2-one)
  517.     SUB    a1, a1, #&10000        ;so arg lies in range [-1,1]
  518.     ADD    a2, a1, a1, LSL #2    ;now apply poly
  519.     ADD    a2, a2, a2, LSL #5
  520.     MOV    a2, a2, ASR #14
  521.     ADD    a2, a2, #&000500
  522.     ADD    a2, a2, #&0000ae
  523.         MOV    ip, a1
  524.     mul16c    a2, ip, a2, a3
  525.     ADD    a2, a2, #&000800
  526.     ADD    a2, a2, #&000088
  527.         MOV    ip, a1
  528.     mul16c    a2, ip, a2, a3
  529.     ADD    a2, a2, #&002000
  530.     ADD    a2, a2, #&00009a
  531.         MOV    ip, a1
  532.     mul16c    a2, ip, a2, a3
  533.     ADD    a2, a2, #&004600
  534.     ADD    a2, a2, #&00009a
  535.         MOV    ip, a1
  536.     mul16c    a2, ip, a2, a3
  537.     ADD    a2, a2, #&030000
  538.     ADD    a2, a2, #&00d900
  539.     ADD    a2, a2, #&0000b4
  540.         MOV    ip, a1
  541.     mul16c    a2, ip, a2, a3
  542.     ADD    a2, a2, #&030000
  543.     ADD    a2, a2, #&00ae00
  544.     ADD    a2, a2, #&000052
  545.     MOV    a1, a2, ASR #4
  546.     TST    a4, #2
  547.     RSBNE    a1, a1, #65536
  548.     TST    a4, #1
  549.     RSBNE    a1, a1, #0
  550.     LDMFD    sp!, {pc}^
  551.  
  552.  
  553.  
  554. ; gauss16
  555. ; a leaf APCS function
  556. ;
  557. ; C prototype:
  558. ; int gauss16(void)
  559. ;
  560. ; returns 65536 * ( pseudo randon variable with approximate distribution N(0,1) )
  561. ; note, the approximate gaussian distribution is achieved via the sum of n pseudo U[0,65535]
  562. ;       random variables & application of the central limit theorem
  563. ;       here we use n=8, giving an actual range of ▒(sqrt24) (*65536).
  564. ;
  565. ;       for general n, recall we need to return:
  566. ;       (sum n U[0,65535] random variables) * 2sqrt(3n)*65536/(65535n)  -  sqrt(3n)*65536
  567. ;
  568.  
  569.     EXPORT    gauss16
  570.  
  571. gansta  DCB     "gauss16", 0
  572.         ALIGN
  573. ganend  DCD     &ff000000 + ganend - gansta
  574.  
  575. gauss16
  576.  
  577.     ADR    a1, gaussseed1
  578.     LDMIA    a1, {a2, a3}
  579.  
  580.     MOVS    a3, a3, LSR #1
  581.         MOVS    a4, a2, RRX
  582.         ADC     a3, a3, a3
  583.         EOR     a4, a4, a2, LSL #12
  584.         EOR     a2, a4, a4, LSR #20
  585.     MOV    ip, a2, LSR #16
  586.     MOV    a4, a2, LSL #16
  587.     ADD    ip, ip, a4, LSR #16
  588.  
  589.     MOVS    a3, a3, LSR #1
  590.         MOVS    a4, a2, RRX
  591.         ADC     a3, a3, a3
  592.         EOR     a4, a4, a2, LSL #12
  593.         EOR     a2, a4, a4, LSR #20
  594.     ADD    ip, ip, a2, LSR #16
  595.     MOV    a4, a2, LSL #16
  596.     ADD    ip, ip, a4, LSR #16
  597.  
  598.     MOVS    a3, a3, LSR #1
  599.         MOVS    a4, a2, RRX
  600.         ADC     a3, a3, a3
  601.         EOR     a4, a4, a2, LSL #12
  602.         EOR     a2, a4, a4, LSR #20
  603.     ADD    ip, ip, a2, LSR #16
  604.     MOV    a4, a2, LSL #16
  605.     ADD    ip, ip, a4, LSR #16
  606.  
  607.     MOVS    a3, a3, LSR #1
  608.         MOVS    a4, a2, RRX
  609.         ADC     a3, a3, a3
  610.         EOR     a4, a4, a2, LSL #12
  611.         EOR     a2, a4, a4, LSR #20
  612.     ADD    ip, ip, a2, LSR #16
  613.     MOV    a4, a2, LSL #16
  614.     ADD    ip, ip, a4, LSR #16        ;sum now in ip
  615.  
  616.     STMIA    a1, {a2, a3}
  617.  
  618.     RSB    a1, ip, ip, ASL #3
  619.     RSB    a2, ip, ip, ASL #2
  620.     ADD    a1, a2, a1, ASL #4
  621.     ADD    a1, a1, ip, ASL #9
  622.     ADD    a2, ip, ip, ASL #4
  623.     ADD    a2, a2, ip, ASL #6
  624.     ADD    a1, a1, a2, LSR #10
  625.     MOV    a1, a1, LSR #9
  626.     SUB    a1, a1, #&4E000
  627.     SUB    a1, a1, #&00620
  628.     SUB    a1, a1, #&00004
  629.  
  630.         MOVS    pc, lr
  631.  
  632. gaussseed1    DCD    -1        ;bits b0-b31 of seed
  633.         DCD    -1        ;bit  b32 of seed (in lsb of word)
  634.  
  635. ; sgauss16
  636. ; a leaf APCS function
  637. ;
  638. ; C prototype:
  639. ; void sgauss16(int seed)
  640. ;
  641. ; sets the seed for gauss16
  642. ;
  643.  
  644.     EXPORT    sgauss16
  645.  
  646. sgnsta  DCB     "sgauss16", 0
  647.         ALIGN
  648. sgnend  DCD     &ff000000 + sgnend - sgnsta
  649.  
  650. sgauss16
  651.  
  652.     ADR    a2, gaussseed1
  653.     MOV    a3, #1
  654.     STMIA    a2, {a1, a3}
  655.         MOVS    pc, lr
  656.  
  657.  
  658.  
  659. ; div_frac16
  660. ; a leaf APCS function
  661. ;
  662. ; C prototype:
  663. ; int div_frac16(int number, int divisor)
  664. ;
  665. ; returns integer part of 65536*number/divisor
  666. ; assumes abs number < 65536 * abs divisor
  667. ; if this needs checking, must be done by caller
  668. ;
  669.  
  670.         EXPORT  div_frac16
  671.  
  672. dfnsta  DCB     "div_frac16", 0
  673.         ALIGN
  674. dfnend  DCD     &ff000000 + dfnend - dfnsta
  675.  
  676. div_frac16
  677.  
  678.         div16   a1, a2, a1, a3, a4
  679.         MOVS    pc, lr
  680.  
  681.  
  682.  
  683. ; mul_frac16
  684. ; a leaf APCS function
  685. ;
  686. ; C prototype:
  687. ; int mul_frac16(int x, int a)
  688. ;
  689. ; returns 32-bit signed integer x*a/65536
  690. ; assumes result fits into 32-bit signed representation
  691. ; note, no other restrictions on a - if can guarantee abs a < 65536, use mul_frac16c instead as is quicker
  692. ;
  693.  
  694.         EXPORT  mul_frac16
  695.  
  696. mfnsta  DCB     "mul_frac16", 0
  697.         ALIGN
  698. mfnend  DCD     &ff000000 + mfnend - mfnsta
  699.  
  700. mul_frac16
  701.  
  702.         mul16   a1, a2, a1, a3, a4, ip
  703.         MOVS    pc, lr
  704.  
  705.  
  706.  
  707. ; mul_frac16c
  708. ; a leaf APCS function
  709. ;
  710. ; C prototype:
  711. ; int mul_frac16c(int x, int a)
  712. ;
  713. ; returns 32-bit signed integer x*a/65536
  714. ; assumes abs a <=65536
  715. ; if it is possible that abs a > 65536, caller must check range & either not call fn or round down to 65536
  716. ;
  717.  
  718.         EXPORT  mul_frac16c
  719.  
  720. mfcnsta DCB     "mul_frac16c", 0
  721.         ALIGN
  722. mfcnend DCD     &ff000000 + mfcnend - mfcnsta
  723.  
  724. mul_frac16c
  725.  
  726.         mul16c  a1, a2, a1, a3
  727.         MOVS    pc, lr
  728.  
  729.  
  730.  
  731. ; sqrt_frac16
  732. ; a leaf APCS function
  733. ;
  734. ; C prototype:
  735. ; int sqrt_frac16(unsigned int x)
  736. ;
  737. ; returns 32-bit integer sqrt(x<<16)
  738. ;
  739.  
  740.         EXPORT  sqrt_frac16
  741.  
  742. sqfnsta DCB     "sqrt_frac16", 0
  743.         ALIGN
  744. sqfnend DCD     &ff000000 + sqfnend - sqfnsta
  745.  
  746. sqrt_frac16
  747.  
  748.         sqrt16  a1, a1, a2, a3, a4, ip
  749.         MOVS    pc, lr
  750.  
  751.  
  752.  
  753.         END
  754.