Unit JFDctFst;

{ This file contains a fast, not so accurate integer implementation of the
  forward DCT (Discrete Cosine Transform).

  A 2-D DCT can be done by 1-D DCT on each row followed by 1-D DCT
  on each column.  Direct algorithms are also available, but they are
  much more complex and seem not to be any faster when reduced to code.

  This implementation is based on Arai, Agui, and Nakajima's algorithm for
  scaled DCT.  Their original paper (Trans. IEICE E-71(11):1095) is in
  Japanese, but the algorithm is described in the Pennebaker & Mitchell
  JPEG textbook (see REFERENCES section in file README).  The following code
  is based directly on figure 4-8 in P&M.
  While an 8-point DCT cannot be done in less than 11 multiplies, it is
  possible to arrange the computation so that many of the multiplies are
  simple scalings of the final outputs.  These multiplies can then be
  folded into the multiplications or divisions by the JPEG quantization
  table entries.  The AA&N method leaves only 5 multiplies and 29 adds
  to be done in the DCT itself.
  The primary disadvantage of this method is that with fixed-point math,
  accuracy is lost due to imprecise representation of the scaled
  quantization values.  The smaller the quantization table entry, the less
  precise the scaled value, so this implementation does worse with high-
  quality-setting files than with low-quality ones. }

{ Original: jfdctfst.c ; Copyright (C) 1994-1996, Thomas G. Lane. }


interface

{$I jconfig.inc}

uses
  jmorecfg,
  jinclude,
  jpeglib,
  jdct;   	{ Private declarations for DCT subsystem }


{ Perform the forward DCT on one block of samples. }

{GLOBAL}
procedure jpeg_fdct_ifast (var data : array of DCTELEM);

implementation

{ This module is specialized to the case DCTSIZE = 8. }

{$ifndef DCTSIZE_IS_8}
  Sorry, this code only copes with 8x8 DCTs. { deliberate syntax err }
{$endif}


{ Scaling decisions are generally the same as in the LL&M algorithm;
  see jfdctint.c for more details.  However, we choose to descale
  (right shift) multiplication products as soon as they are formed,
  rather than carrying additional fractional bits into subsequent additions.
  This compromises accuracy slightly, but it lets us save a few shifts.
  More importantly, 16-bit arithmetic is then adequate (for 8-bit samples)
  everywhere except in the multiplications proper; this saves a good deal
  of work on 16-bit-int machines.

  Again to save a few shifts, the intermediate results between pass 1 and
  pass 2 are not upscaled, but are represented only to integral precision.

  A final compromise is to represent the multiplicative constants to only
  8 fractional bits, rather than 13.  This saves some shifting work on some
  machines, and may also reduce the cost of multiplication (since there
  are fewer one-bits in the constants). }

const
  CONST_BITS = 8;
const
  CONST_SCALE = (INT32(1) shl CONST_BITS);


const
  FIX_0_382683433 = INT32(Round(CONST_SCALE * 0.382683433)); {98}
  FIX_0_541196100 = INT32(Round(CONST_SCALE * 0.541196100)); {139}
  FIX_0_707106781 = INT32(Round(CONST_SCALE * 0.707106781)); {181}
  FIX_1_306562965 = INT32(Round(CONST_SCALE * 1.306562965)); {334}

{ Descale and correctly round an INT32 value that's scaled by N bits.
  We assume RIGHT_SHIFT rounds towards minus infinity, so adding
  the fudge factor is correct for either sign of X. }

function DESCALE(x : INT32; n : int) : INT32;
var
  shift_temp : INT32;
begin
{ We can gain a little more speed, with a further compromise in accuracy,
  by omitting the addition in a descaling shift.  This yields an incorrectly
  rounded result half the time... }
{$ifndef USE_ACCURATE_ROUNDING}
  shift_temp := x;
{$else}
  shift_temp := x + (INT32(1) shl (n-1));
{$endif}

{$ifdef RIGHT_SHIFT_IS_UNSIGNED}
  if shift_temp < 0 then
    Descale :=  (shift_temp shr n) or ((not INT32(0)) shl (32-n))
  else
{$endif}
    Descale :=  (shift_temp shr n);
end;

{ Multiply a DCTELEM variable by an INT32 constant, and immediately
  descale to yield a DCTELEM result. }


   function MULTIPLY(X : DCTELEM; Y: INT32): DCTELEM;
   begin
     Multiply := DeScale((X) * (Y), CONST_BITS);
   end;


{ Perform the forward DCT on one block of samples. }

{GLOBAL}
procedure jpeg_fdct_ifast (var data : array of DCTELEM);
type
  PWorkspace = ^TWorkspace;
  TWorkspace = array [0..DCTSIZE2-1] of DCTELEM;
var
  tmp0, tmp1, tmp2, tmp3, tmp4, tmp5, tmp6, tmp7 : DCTELEM;
  tmp10, tmp11, tmp12, tmp13 : DCTELEM;
  z1, z2, z3, z4, z5, z11, z13 : DCTELEM;
  dataptr :  PWorkspace;
  ctr : int;
  {SHIFT_TEMPS}
begin
  { Pass 1: process rows. }

  dataptr := PWorkspace(@data);
  for ctr := DCTSIZE-1 downto 0 do
  begin
    tmp0 := dataptr^[0] + dataptr^[7];
    tmp7 := dataptr^[0] - dataptr^[7];
    tmp1 := dataptr^[1] + dataptr^[6];
    tmp6 := dataptr^[1] - dataptr^[6];
    tmp2 := dataptr^[2] + dataptr^[5];
    tmp5 := dataptr^[2] - dataptr^[5];
    tmp3 := dataptr^[3] + dataptr^[4];
    tmp4 := dataptr^[3] - dataptr^[4];

    { Even part }

    tmp10 := tmp0 + tmp3;	{ phase 2 }
    tmp13 := tmp0 - tmp3;
    tmp11 := tmp1 + tmp2;
    tmp12 := tmp1 - tmp2;

    dataptr^[0] := tmp10 + tmp11; { phase 3 }
    dataptr^[4] := tmp10 - tmp11;

    z1 := MULTIPLY(tmp12 + tmp13, FIX_0_707106781); { c4 }
    dataptr^[2] := tmp13 + z1;	{ phase 5 }
    dataptr^[6] := tmp13 - z1;

    { Odd part }

    tmp10 := tmp4 + tmp5;	{ phase 2 }
    tmp11 := tmp5 + tmp6;
    tmp12 := tmp6 + tmp7;

    { The rotator is modified from fig 4-8 to avoid extra negations. }
    z5 := MULTIPLY(tmp10 - tmp12, FIX_0_382683433); { c6 }
    z2 := MULTIPLY(tmp10, FIX_0_541196100) + z5; { c2-c6 }
    z4 := MULTIPLY(tmp12, FIX_1_306562965) + z5; { c2+c6 }
    z3 := MULTIPLY(tmp11, FIX_0_707106781); { c4 }

    z11 := tmp7 + z3;		{ phase 5 }
    z13 := tmp7 - z3;

    dataptr^[5] := z13 + z2;	{ phase 6 }
    dataptr^[3] := z13 - z2;
    dataptr^[1] := z11 + z4;
    dataptr^[7] := z11 - z4;

    Inc(DCTELEMPTR(dataptr), DCTSIZE);	{ advance pointer to next row }
  end;

  { Pass 2: process columns. }

  dataptr := PWorkspace(@data);
  for ctr := DCTSIZE-1 downto 0 do
  begin
    tmp0 := dataptr^[DCTSIZE*0] + dataptr^[DCTSIZE*7];
    tmp7 := dataptr^[DCTSIZE*0] - dataptr^[DCTSIZE*7];
    tmp1 := dataptr^[DCTSIZE*1] + dataptr^[DCTSIZE*6];
    tmp6 := dataptr^[DCTSIZE*1] - dataptr^[DCTSIZE*6];
    tmp2 := dataptr^[DCTSIZE*2] + dataptr^[DCTSIZE*5];
    tmp5 := dataptr^[DCTSIZE*2] - dataptr^[DCTSIZE*5];
    tmp3 := dataptr^[DCTSIZE*3] + dataptr^[DCTSIZE*4];
    tmp4 := dataptr^[DCTSIZE*3] - dataptr^[DCTSIZE*4];

    { Even part }

    tmp10 := tmp0 + tmp3;	{ phase 2 }
    tmp13 := tmp0 - tmp3;
    tmp11 := tmp1 + tmp2;
    tmp12 := tmp1 - tmp2;

    dataptr^[DCTSIZE*0] := tmp10 + tmp11; { phase 3 }
    dataptr^[DCTSIZE*4] := tmp10 - tmp11;

    z1 := MULTIPLY(tmp12 + tmp13, FIX_0_707106781); { c4 }
    dataptr^[DCTSIZE*2] := tmp13 + z1; { phase 5 }
    dataptr^[DCTSIZE*6] := tmp13 - z1;

    { Odd part }

    tmp10 := tmp4 + tmp5;	{ phase 2 }
    tmp11 := tmp5 + tmp6;
    tmp12 := tmp6 + tmp7;

    { The rotator is modified from fig 4-8 to avoid extra negations. }
    z5 := MULTIPLY(tmp10 - tmp12, FIX_0_382683433); { c6 }
    z2 := MULTIPLY(tmp10, FIX_0_541196100) + z5; { c2-c6 }
    z4 := MULTIPLY(tmp12, FIX_1_306562965) + z5; { c2+c6 }
    z3 := MULTIPLY(tmp11, FIX_0_707106781); { c4 }

    z11 := tmp7 + z3;		{ phase 5 }
    z13 := tmp7 - z3;

    dataptr^[DCTSIZE*5] := z13 + z2; { phase 6 }
    dataptr^[DCTSIZE*3] := z13 - z2;
    dataptr^[DCTSIZE*1] := z11 + z4;
    dataptr^[DCTSIZE*7] := z11 - z4;

    Inc(DCTELEMPTR(dataptr));	{ advance pointer to next column }
  end;
end;

end.