home *** CD-ROM | disk | FTP | other *** search
/ AmigActive 13 / AACD13.ISO / AACD / Sound / LAME / src / i386 / fftsse.nas < prev    next >
Text File  |  2000-01-30  |  16KB  |  539 lines
  1. ; back port from GOGO-no coda 2.24b by Takehiro TOMINAGA
  2.  
  3. ; GOGO-no-coda
  4. ;    Copyright (C) 1999 shigeo
  5. ;    special thanks to Keiichi SAKAI
  6.  
  7. %include "nasm.h"
  8.  
  9.     globaldef fht_SSE
  10.     globaldef fft_side_SSE
  11.     externdef costab_fft
  12.     externdef sintab_fft
  13.  
  14.     segment_data
  15.     align 16
  16. Q_MMPP    dd    0x0,0x0,0x80000000,0x80000000
  17. Q_MPMP    dd    0x0,0x80000000,0x0,0x80000000
  18. Q_002    dd    0.02236068, 0.02236068, 0.02236068, 0.02236068
  19. D_SQRT2    dd    1.414213562,1.414213562
  20. S_025    dd    0.25
  21. S_05    DD    0.5
  22. S_00005    DD    0.0005
  23.  
  24.     segment_code
  25. ;------------------------------------------------------------------------
  26. ;    by K. SAKAI
  27. ;    99/08/18    PIII 23k[clk]
  28. ;    99/08/19    Ì¿Îá½ç½øÆþ¤ì´¹¤¨ PIII 22k[clk]
  29. ;    99/08/20    bit reversal ¤òµì¸á¸å¤«¤é°Ü¿¢¤·¤¿ PIII 17k[clk]
  30. ;    99/08/23    °ìÉô unroll PIII 14k[clk]
  31. ;    99/11/12    clean up
  32. ;
  33. ;void fht_SSE(float *fz, int n);
  34.     align 16
  35. fht_SSE:
  36.     push    ebx
  37.     push    esi
  38.     push    edi
  39.     push    ebp
  40. %assign _P 4*4
  41.  
  42.     ;2¤ÄÌܤΥ롼¥×
  43.     mov    eax,[esp+_P+4]    ;eax=fz
  44.     mov    ebp,[esp+_P+8]    ;=n
  45.     shl    ebp,2
  46.     add    ebp,eax        ; fn  = fz + n, ¤³¤Î´Ø¿ô½ªÎ»¤Þ¤ÇÉÔÊÑ
  47.  
  48.     xor    ecx,ecx        ; ecx=k=0
  49.     xor    eax,eax
  50.     mov    al,4        ; =k1=1*(sizeof float)    // 4, 16, 64, 256,...
  51.     xor    edx,edx
  52.     mov    dl,12        ; =k3=3*k1
  53.     jmp    short .lp2
  54.  
  55.     align    16
  56. .lp2:                ; do{
  57.     add    cl,2        ; k  += 2;
  58.     shl    eax,2
  59.     shl    edx,2
  60.  
  61.     mov    esi,[esp+_P+4]    ;esi=fi=fz
  62.     mov    edi,eax
  63.     shr    edi,1
  64.     add    edi,esi        ; edi=gi=fi+ki/2
  65.  
  66. ; ¤¿¤«¤À¤«2ÊÂÎó¤·¤«´üÂԤǤ­¤Ê¤¤Éôʬ¤ÏFPU¤Î¤Û¤¦¤¬Â®¤¤¡£
  67.     movss    xmm7,[D_SQRT2]
  68.     jmp    short .lp20
  69.  
  70.     align    16
  71. .lp20:                ; do{
  72. ;                       f0     = fi[0 ] + fi[k1];
  73. ;                       f2     = fi[k2] + fi[k3];
  74. ;                       f1     = fi[0 ] - fi[k1];
  75. ;                       f3     = fi[k2] - fi[k3];
  76. ;                       fi[0 ] = f0     + f2;
  77. ;                       fi[k1] = f1     + f3;
  78. ;                       fi[k2] = f0     - f2;
  79. ;                       fi[k3] = f1     - f3;
  80.     fld    dword [esi]
  81.     fadd    dword [esi+eax]
  82.     fld    dword [esi+eax*2]
  83.     fadd    dword [esi+edx]
  84.  
  85.     fld    dword [esi]
  86.     fsub    dword [esi+eax]
  87.     fld    dword [esi+eax*2]
  88.     fsub    dword [esi+edx]
  89.  
  90.     fld    st1
  91.     fadd    st0,st1
  92.     fstp    dword [esi+eax]
  93.     fsubp    st1,st0
  94.     fstp    dword [esi+edx]
  95.  
  96.     fld    st1
  97.     fadd    st0,st1
  98.     fstp    dword [esi]
  99.     fsubp    st1,st0
  100.     fstp    dword [esi+eax*2]
  101.  
  102.     lea    esi,[esi + eax*4]    ; = fi += (k1 * 4);
  103. ;    add    esi,eax
  104. ;    add    esi,edx
  105. ;                       g0     = gi[0 ] + gi[k1];
  106. ;                       g2     = SQRT2  * gi[k2];
  107. ;                       g1     = gi[0 ] - gi[k1];
  108. ;                       g3     = SQRT2  * gi[k3];
  109. ;                       gi[0 ] = g0     + g2;
  110. ;                       gi[k2] = g0     - g2;
  111. ;                       gi[k1] = g1     + g3;
  112. ;                       gi[k3] = g1     - g3;
  113.     fld    dword [edi]
  114.     fadd    dword [edi+eax]
  115.     fld    dword [D_SQRT2]
  116.     fmul    dword [edi+eax*2]
  117.  
  118.     fld    dword [edi]
  119.     fsub    dword [edi+eax]
  120.     fld    dword [D_SQRT2]
  121.     fmul    dword [edi+edx]
  122.  
  123.     fld    st1
  124.     fadd    st0,st1
  125.     fstp    dword [edi+eax]
  126.     fsubp    st1,st0
  127.     fstp    dword [edi+edx]
  128.  
  129.     fld    st1
  130.     fadd    st0,st1
  131.     fstp    dword [edi]
  132.     fsubp    st1,st0
  133.     fstp    dword [edi+eax*2]
  134.  
  135.     lea    edi,[edi + eax*4]    ; = gi += (k1 * 4);
  136.     cmp    esi,ebp
  137.     jl    near .lp20        ; while (fi<fn);
  138.  
  139. ;               i = 1; //for (i=1;i<kx;i++){
  140. ;                       c1 = 1.0*t_c - 0.0*t_s;
  141. ;                       s1 = 0.0*t_c + 1.0*t_s;
  142.     movss    xmm6,[costab_fft + ecx*4]
  143.     movss    xmm1,[sintab_fft + ecx*4]
  144.     shufps    xmm6,xmm1,0x00    ; = {s1, s1, c1, c1}
  145.     shufps    xmm6,xmm6,0x28    ; = {+c1, +s1, +s1, +c1}
  146. ;                       c2 = c1*c1 - s1*s1;
  147. ;                       s2 = c1*s1 + s1*c1;
  148.     movaps    xmm4,xmm6
  149.     movaps    xmm7,xmm6
  150.     unpcklps    xmm4,xmm4    ; = {s1, s1, c1, c1}
  151.     shufps    xmm7,xmm7,0x41
  152.     mulps    xmm4,xmm6    ; = {s1*c1, s1*s1, c1*s1, c1*c1}
  153.     xorps    xmm7,[Q_MMPP]    ; = {-s1, -c1, +c1, +s1}
  154.     movhlps    xmm3,xmm4
  155.     xorps    xmm3,[Q_MPMP]
  156.     subps    xmm4,xmm3    ; = {--, --, s2, c2}
  157.     movlhps    xmm4,xmm4    ; = {+s2, +c2, +s2, +c2}
  158.     movaps    xmm5,xmm4
  159.     shufps    xmm5,xmm5,0x11
  160.     xorps    xmm5,[Q_MPMP]    ; = {-c2, +s2, -c2, +s2}
  161.     mov    esi,[esp+_P+4]    ; = fz
  162.     lea    edi,[esi + eax - 4]    ; edi = gi = fz +k1-i
  163.     add    esi,4        ; esi = fi = fz + i
  164.     jmp    short .lp21
  165.  
  166.     align    16
  167. .lp21:                ; do{
  168. ;                               a       = c2*fi[k1] + s2*gi[k1];
  169. ;                               b       = s2*fi[k1] - c2*gi[k1];
  170. ;                               c       = c2*fi[k3] + s2*gi[k3];
  171. ;                               d       = s2*fi[k3] - c2*gi[k3];
  172. ;                               f0      = fi[0 ]        + a;
  173. ;                               g0      = gi[0 ]        + b;
  174. ;                               f2      = fi[k1 * 2]    + c;
  175. ;                               g2      = gi[k1 * 2]    + d;
  176. ;                               f1      = fi[0 ]        - a;
  177. ;                               g1      = gi[0 ]        - b;
  178. ;                               f3      = fi[k1 * 2]    - c;
  179. ;                               g3      = gi[k1 * 2]    - d;
  180.  
  181.     movss    xmm0,[esi + eax]    ; = fi[k1]
  182.     movss    xmm2,[esi + edx]    ; = fi[k3]
  183.     shufps    xmm0,xmm2,0x00    ; = {fi[k3], fi[k3], fi[k1], fi[k1]}
  184.     movss    xmm1,[edi + eax]    ; = fi[k1]
  185.     movss    xmm3,[edi + edx]    ; = fi[k3]
  186.     shufps    xmm1,xmm3,0x00    ; = {gi[k3], gi[k3], gi[k1], gi[k1]}
  187.     movss    xmm2,[esi]        ; = fi[0]
  188.     mulps    xmm0,xmm4        ; *= {+s2, +c2, +s2, +c2}
  189.     movss    xmm3,[esi + eax*2]    ; = fi[k2]
  190.     unpcklps    xmm2,xmm3    ; = {--, --, fi[k2], fi[0]}
  191.     mulps    xmm1,xmm5        ; *= {-c2, +s2, -c2, +s2}
  192.     movss    xmm3,[edi + eax*2]    ; = gi[k2]
  193.     addps    xmm0,xmm1        ; = {d, c, b, a}
  194.     movss    xmm1,[edi]        ; = gi[0]
  195.     unpcklps    xmm1,xmm3    ; = {--,  --, gi[k2], gi[0]}
  196.     unpcklps    xmm2,xmm1    ; = {gi[k2], fi[k2], gi[0], fi[0]}
  197.     movaps    xmm1,xmm2
  198.     addps    xmm1,xmm0    ; = {g2, f2, g0, f0}
  199.     subps    xmm2,xmm0    ; = {g3, f3, g1, f1}
  200.  
  201. ;                               a       = c1*f2     + s1*g3;
  202. ;                               c       = s1*g2     + c1*f3;
  203. ;                               b       = s1*f2     - c1*g3;
  204. ;                               d       = c1*g2     - s1*f3;
  205. ;                               fi[0 ]  = f0        + a;
  206. ;                               gi[0 ]  = g0        + c;
  207. ;                               gi[k1]  = g1        + b;
  208. ;                               fi[k1]  = f1        + d;
  209. ;                               fi[k1 * 2]  = f0    - a;
  210. ;                               gi[k1 * 2]  = g0    - c;
  211. ;                               gi[k3]      = g1    - b;
  212. ;                               fi[k3]      = f1    - d;
  213.     movaps    xmm3,xmm1
  214.     movhlps    xmm1,xmm1    ; = {g2, f2, g2, f2}
  215.     shufps    xmm3,xmm2,0x14    ; = {f1, g1, g0, f0}
  216.     mulps    xmm1,xmm6    ; *= {+c1, +s1, +s1, +c1}
  217.     shufps    xmm2,xmm2,0xBB    ; = {f3, g3, f3, g3}
  218.     mulps    xmm2,xmm7    ; *= {-s1, -c1, +c1, +s1}
  219.     addps    xmm1,xmm2    ; = {d, b, c, a}
  220.     movaps    xmm2,xmm3
  221.     addps    xmm3,xmm1    ; = {fi[k1], gi[k1], gi[0], fi[0]}
  222.     subps    xmm2,xmm1    ; = {fi[k3], gi[k3], gi[k1*2], fi[k1*2]}
  223.     movhlps    xmm0,xmm3
  224.     movss    [esi],xmm3
  225.     shufps    xmm3,xmm3,0x55
  226.     movss    [edi+eax],xmm0
  227.     shufps    xmm0,xmm0,0x55
  228.     movss    [edi],xmm3
  229.     movss    [esi+eax],xmm0
  230.     movhlps    xmm0,xmm2
  231.     movss    [esi+eax*2],xmm2
  232.     shufps    xmm2,xmm2,0x55
  233.     movss    [edi+edx],xmm0
  234.     shufps    xmm0,xmm0,0x55
  235.     movss    [edi+eax*2],xmm2
  236.     lea    edi,[edi + eax*4] ; gi += (k1 * 4);
  237.     movss    [esi+edx],xmm0
  238.     lea    esi,[esi + eax*4] ; fi += (k1 * 4);
  239.     cmp    esi,ebp
  240.     jl    near .lp21        ; while (fi<fn);
  241. ; unrollÁ°¤Îdo loop¤Ï43+4Ì¿Îá
  242.  
  243. ; ºÇÆâ¼þ¤Ç¤Ï¤Ê¤¤for¥ë¡¼¥×¤òunrolling¤·¤¿
  244. ; kx=   2,   8,  32,  128
  245. ; k4=  16,  64, 256, 1024
  246. ;       0, 6/2,30/2,126/2
  247. ; at here
  248. ;    xmm6 = {--, --, s1, c1}
  249. ;               c3 = c1; s3 = s1;
  250.     xor    ebx,ebx
  251.     mov    bl,4*4        ; = i = 4
  252.     cmp    ebx,eax        ; i < k1
  253.     jnl    near .F22
  254.  
  255.     shufps    xmm6,xmm6,0x14    ; = {c1, s1, s1, c1}
  256.     jmp    short .F220
  257.  
  258.     align    16
  259. ;               for (i=4;i<k1;i+=4){ // for (i=2;i<k1/2;i+=2){
  260. .lp22:
  261.     shufps    xmm6,xmm6,0x69    ; xmm6 = {c3, s3, s3, c3}
  262. .F220:
  263. ; at here, xmm6 is {c3, s3, s3, c3}
  264. ;                       c1 = c3*t_c - s3*t_s;
  265. ;                       s1 = c3*t_s + s3*t_c;
  266.     movss    xmm0,[costab_fft + ecx*4]
  267.     movss    xmm1,[sintab_fft + ecx*4]
  268.     shufps    xmm0,xmm1,0x00    ; = {t_s, t_s, t_c, t_c}
  269.     mulps    xmm6,xmm0
  270.     movhlps    xmm4,xmm6
  271.     xorps    xmm4,[Q_MPMP]
  272.     subps    xmm6,xmm4    ; = {--, --, s1, c1}
  273.  
  274. ;                       c3 = c1*t_c - s1*t_s;
  275. ;                       s3 = s1*t_c + c1*t_s;
  276.     shufps    xmm6,xmm6,0x14    ; = {c1, s1, s1, c1}
  277.     mulps    xmm0,xmm6
  278.     movhlps    xmm3,xmm0
  279.     xorps    xmm3,[Q_MPMP]
  280.     subps    xmm0,xmm3    ; = {--, --, s3, c3}
  281.  
  282.     unpcklps    xmm6,xmm0    ; xmm6 = {s3, s1, c3, c1}
  283.     shufps    xmm6,xmm6,0xB4    ; xmm6 = {s1, s3, c3, c1}
  284.  
  285. ;                       c2 = c1*c1 - s1*s1;
  286. ;                       c4 = c3*c3 - s3*s3;
  287. ;                       s4 = s3*c3 + s3*c3;
  288. ;                       s2 = s1*c1 + s1*c1;
  289.     movaps    xmm7,xmm6
  290.     movaps    xmm4,xmm6
  291.     shufps    xmm7,xmm7,0x14
  292.     shufps    xmm4,xmm4,0xEB
  293.     xorps    xmm4,[Q_MMPP]    ; = {-c3,-c1, s3, s1}
  294.     mulps    xmm7,xmm6
  295.     mulps    xmm4,xmm6
  296.     shufps    xmm4,xmm4,0x1B
  297.     addps    xmm7,xmm4    ; xmm7 = {s2, s4, c4, c2}
  298.  
  299. ;                       fi = fz +i;
  300. ;                       gi = fz +k1-i;
  301.     mov    edi,[esp+_P+4]    ; = fz
  302.     mov    esi,ebx
  303.     shr    esi,1
  304.     sub    edi,esi        ; edi = fz - i/2
  305.     lea    esi,[edi + ebx]    ; esi = fi = fz +i/2
  306.     add    edi,eax        ; edi = gi = fz +k1-i/2
  307.     sub    edi,4
  308. ;                       do{
  309. .lp220:
  310. ; unroll¸å¤Îdo loop¤Ï51+4Ì¿Îá
  311. ;                               a       = c2*fi[k1  ] + s2*gi[k1  ];
  312. ;                               e       = c4*fi[k1+1] + s4*gi[k1-1];
  313. ;                               f       = s4*fi[k1+1] - c4*gi[k1-1];
  314. ;                               b       = s2*fi[k1  ] - c2*gi[k1  ];
  315. ;                               c       = c2*fi[k3  ] + s2*gi[k3  ];
  316. ;                               g       = c4*fi[k3+1] + s4*gi[k3-1];
  317. ;                               h       = s4*fi[k3+1] - c4*gi[k3-1];
  318. ;                               d       = s2*fi[k3  ] - c2*gi[k3  ];
  319.  
  320.     movaps    xmm4,xmm7    ; xmm7 = {s2, s4, c4, c2}
  321.     shufps    xmm4,xmm4,0x1B
  322.     xorps    xmm4,[Q_MMPP]
  323.     movlps    xmm0,[esi+eax]
  324.     movlps    xmm1,[edi+eax]
  325.     movlps    xmm2,[esi+edx]
  326.     movlps    xmm3,[edi+edx]
  327.     shufps    xmm0,xmm0,0x14
  328.     shufps    xmm1,xmm1,0x41
  329.     shufps    xmm2,xmm2,0x14
  330.     shufps    xmm3,xmm3,0x41
  331.     mulps    xmm0,xmm7
  332.     mulps    xmm1,xmm4
  333.     mulps    xmm2,xmm7
  334.     mulps    xmm3,xmm4
  335.     addps    xmm0,xmm1    ; xmm0 = {b, f, e, a}
  336.     addps    xmm2,xmm3    ; xmm2 = {d, h, g, c}
  337. ;17
  338.  
  339. ;                               f0      = fi[0   ]    + a;
  340. ;                               f4      = fi[0 +1]    + e;
  341. ;                               g4      = gi[0 -1]    + f;
  342. ;                               g0      = gi[0   ]    + b;
  343. ;                               f1      = fi[0   ]    - a;
  344. ;                               f5      = fi[0 +1]    - e;
  345. ;                               g5      = gi[0 -1]    - f;
  346. ;                               g1      = gi[0   ]    - b;
  347. ;                               f2      = fi[k2  ]    + c;
  348. ;                               f6      = fi[k2+1]    + g;
  349. ;                               g6      = gi[k2-1]    + h;
  350. ;                               g2      = gi[k2  ]    + d;
  351. ;                               f3      = fi[k2  ]    - c;
  352. ;                               f7      = fi[k2+1]    - g;
  353. ;                               g7      = gi[k2-1]    - h;
  354. ;                               g3      = gi[k2  ]    - d;
  355.     movlps    xmm4,[esi      ]
  356.     movhps    xmm4,[edi      ]
  357.     movaps    xmm1,xmm4
  358.     subps    xmm1,xmm0    ; xmm1 = {g1, g5, f5, f1}
  359.     movlps    xmm5,[esi+eax*2]
  360.     movhps    xmm5,[edi+eax*2]
  361.     movaps    xmm3,xmm5
  362.     subps    xmm3,xmm2    ; xmm3 = {g3, g7, f7, f3}
  363.     addps    xmm0,xmm4    ; xmm0 = {g0, g4, f4, f0}
  364.     addps    xmm2,xmm5    ; xmm2 = {g2, g6, f6, f2}
  365. ;10
  366.  
  367. ;                               a       = c1*f2     + s1*g3;    ½ç*½ç + µÕ*µÕ
  368. ;                               e       = c3*f6     + s3*g7;
  369. ;                               g       = s3*g6     + c3*f7;
  370. ;                               c       = s1*g2     + c1*f3;
  371. ;                               d       = c1*g2     - s1*f3;    ½ç*µÕ - µÕ*½ç
  372. ;                               h       = c3*g6     - s3*f7;
  373. ;                               f       = s3*f6     - c3*g7;
  374. ;                               b       = s1*f2     - c1*g3;
  375.  
  376.     movaps    xmm5,xmm6    ; xmm6 = {s1, s3, c3, c1}
  377.     shufps    xmm5,xmm5,0x1B    ; = {c1, c3, s3, s1}
  378.     movaps    xmm4,xmm2
  379.     mulps    xmm4,xmm6
  380.     shufps    xmm2,xmm2,0x1B    ; xmm2 = {f2, f6, g6, g2}
  381.     mulps    xmm2,xmm6
  382.     mulps    xmm5,xmm3
  383.     mulps    xmm3,xmm6
  384.     shufps    xmm3,xmm3,0x1B
  385.     addps    xmm4,xmm3    ; = {c, g, e, a}
  386.     subps    xmm2,xmm5    ; = {b, f, h, d}
  387. ;10
  388.  
  389. ;                               fi[0   ]  = f0        + a;
  390. ;                               fi[0 +1]  = f4        + e;
  391. ;                               gi[0 -1]  = g4        + g;
  392. ;                               gi[0   ]  = g0        + c;
  393. ;                               fi[k2  ]  = f0        - a;
  394. ;                               fi[k2+1]  = f4        - e;
  395. ;                               gi[k2-1]  = g4        - g;
  396. ;                               gi[k2  ]  = g0        - c;
  397. ;                               fi[k1  ]  = f1        + d;
  398. ;                               fi[k1+1]  = f5        + h;
  399. ;                               gi[k1-1]  = g5        + f;
  400. ;                               gi[k1  ]  = g1        + b;
  401. ;                               fi[k3  ]  = f1        - d;
  402. ;                               fi[k3+1]  = f5        - h;
  403. ;                               gi[k3-1]  = g5        - f;
  404. ;                               gi[k3  ]  = g1        - b;
  405.     movaps    xmm3,xmm0
  406.     subps    xmm0,xmm4
  407.     movlps    [esi+eax*2],xmm0
  408.     movhps    [edi+eax*2],xmm0
  409.     addps    xmm4,xmm3
  410.     movlps    [esi      ],xmm4
  411.     movhps    [edi      ],xmm4
  412.  
  413.     movaps    xmm5,xmm1
  414.     subps    xmm1,xmm2
  415.     movlps    [esi+edx  ],xmm1
  416.     movhps    [edi+edx  ],xmm1
  417.     addps    xmm2,xmm5
  418.     movlps    [esi+eax  ],xmm2
  419.     movhps    [edi+eax  ],xmm2
  420. ; 14
  421. ;                               gi     += k4;
  422. ;                               fi     += k4;
  423.     lea    edi,[edi + eax*4] ; gi += (k1 * 4);
  424.     lea    esi,[esi + eax*4] ; fi += (k1 * 4);
  425.     cmp    esi,ebp
  426.     jl    near .lp220        ; while (fi<fn);
  427. ;                       } while (fi<fn);
  428.  
  429.     add    ebx,4*4        ; i+= 4
  430.     cmp    ebx,eax        ; i < k1
  431.     jl    near .lp22
  432. ;               }
  433. .F22:
  434.  
  435.     cmp    eax,[esp+_P+8]    ; while ((k1 * 4)<n);
  436.     jl    near .lp2
  437.  
  438.     pop    ebp
  439.     pop    edi
  440.     pop    esi
  441.     pop    ebx
  442.     ret
  443.  
  444. ;------------------------------------------------------------------------
  445. ;    99/11/12    Initial version for SSE by K. SAKAI, 4300clk@P3
  446. ; This routine is very slow when wsamp_r_int is not aligned to 16byte boundary.
  447. ;
  448. ;void fft_side_SSE( float in[2][1024], int s, float *ret)
  449. ;        energy = (in[0][512] - in[1][512])^2;
  450. ;        energy = (in[0][1024-s] - in[1][1024-s])^2;
  451. ;        for (i=s,j=1024-s;i<512;i++,j--){
  452. ;                a = in[0][i] - in[1][i];
  453. ;                energy += a*a;
  454. ;                b = in[0][j-1] - in[1][j-1];
  455. ;                energy += b*b;
  456. ;        }
  457. ;        *ret = energy * 0.25;
  458.  
  459.     align    16
  460. fft_side_SSE:
  461.     mov    ecx,[esp+8]    ; = i = s
  462.     mov    edx,1024
  463.     sub    edx,ecx        ; = j
  464.     mov    eax,[esp+4]    ; = in
  465.     movss    xmm7,[eax+1024*0*4+512*4]
  466.     movss    xmm1,[eax+1024*1*4+512*4]
  467.     subss    xmm7,xmm1
  468.     mulss    xmm7,xmm7
  469.     movss    xmm2,[eax+1024*0*4+edx*4]
  470.     movss    xmm3,[eax+1024*1*4+edx*4]
  471.     subss    xmm2,xmm3
  472.     mulss    xmm2,xmm2
  473.     addss    xmm7,xmm2
  474.  
  475.     test    cl,1
  476.     jz    .even
  477.  
  478. .odd:    dec    edx
  479.     movss    xmm0,[eax+1024*0*4+ecx*4]
  480.     movss    xmm1,[eax+1024*1*4+ecx*4]
  481.     inc    ecx
  482.     movss    xmm2,[eax+1024*0*4+edx*4]
  483.     movss    xmm3,[eax+1024*1*4+edx*4]
  484.     cmp    ecx,edx
  485.     subss    xmm0,xmm1
  486.     subss    xmm2,xmm3
  487.     mulss    xmm0,xmm0
  488.     mulss    xmm2,xmm2
  489.     addss    xmm7,xmm0
  490.     addss    xmm7,xmm2
  491.     je    near .exit1
  492.  
  493. .even:    test    cl,2
  494.     jz    .f0
  495.     sub    edx,2
  496.     movlps    xmm0,[eax+1024*0*4+ecx*4]
  497.     movlps    xmm1,[eax+1024*1*4+ecx*4]
  498.     add    ecx,2
  499.     movhps    xmm0,[eax+1024*0*4+edx*4]
  500.     movhps    xmm1,[eax+1024*1*4+edx*4]
  501.     cmp    ecx,edx
  502.     subps    xmm0,xmm1
  503.     mulps    xmm0,xmm0
  504.     addps    xmm7,xmm0
  505.     je    .exit4
  506.     jmp    short .f0
  507.  
  508.     align    16
  509. .f0:
  510. .lp0:
  511.     sub    edx,4
  512.     movaps    xmm0,[eax+1024*0*4+ecx*4]
  513.     movaps    xmm1,[eax+1024*1*4+ecx*4]
  514.     add    ecx,4
  515.     subps    xmm0,xmm1
  516.     mulps    xmm0,xmm0
  517.     addps    xmm7,xmm0
  518.     movaps    xmm2,[eax+1024*0*4+edx*4]
  519.     movaps    xmm3,[eax+1024*1*4+edx*4]
  520.     cmp    ecx,edx
  521.     subps    xmm2,xmm3
  522.     mulps    xmm2,xmm2
  523.     addps    xmm7,xmm2
  524.     jne    .lp0
  525.  
  526. .exit4:    movhlps    xmm6,xmm7
  527.     addps    xmm7,xmm6
  528.     movaps    xmm6,xmm7
  529.     shufps    xmm6,xmm6,01010101B
  530.     addss    xmm7,xmm6
  531.  
  532. .exit1:    mulss    xmm7,[S_025]
  533.     mov    eax,[esp+12]
  534.     movss    [eax],xmm7
  535.     ret
  536.  
  537.  
  538.     end
  539.