Форум » [x]Harbour » И снова Round » Ответить
И снова Round
NickSam: Нужно округление с большой точностью, стандартный Round как известно глючит. Кто может поделиться исходниками какого нибудь нормального Round ? Поиск по форуму дает вот эту тему, но там уже все ссылки нерабочие: [url=https://clipper.borda.ru/?1-0-0-00000430-000-0-0-1320071469]Кто может поделиться исходниками RoundTo(N,Dec,[Round5toUp] от suv7[/url]
Ответов - 10
Dima: [pre2] Function P_Round(Number,Decimals) static Const := 0.0000000000000002220446049250313 Decimals:=if(Decimals==NIL,2,Decimals) return round((mantissa(Number)+if(Number>0,1,-1)*Const)*2^exponent(Number),Decimals) [/pre2] А еще есть Ft_round
Dima: о нашел , счас
Dima: [pre2] #include "suv.ch" // FloatDigits -количество значащих позиций, принимающихся во внимание // при операциях сравнения, округления и т.п. Значащие позиции расходуются // не только под дробную часть числа, но и под целую // Устанавливать меньше 12 не стоит // Устанавливать больше 12 нет смысла, т.к. при обычных операциях может накапливаться ошибка уже за 12м знаком (см пример в комментариях с 0.5035) // Устанавливать больше 15-16 нет тем более никакого смысла, т.к. это больше макс знач позиций в double // 12 позиций обеспечивают нормальное округление сотен миллионов до копеек // при этом точность будет составлять +-0.0005 (для сотен миллионов) // 123 456 789.012 static FloatDigits:=12 static Round5DefaultMode:=.T. // .T.=5 вверх .F.=округление к четному setget func FloatDigits setget func Round5toUp for Round5DefaultMode *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func SignifPower(N) // значащая степень local i,IntDigits,SignificantPower:=0,Sign IntDigits := len(nTrim(int(N))) // количество позиций в целой части числа if IntDigits < FloatDigits // если целая часть занимает меньше позиций, чем задано в FloatDigits // то некоторое кол-во дробных позиций являются значимыми SignificantPower := -(FloatDigits - IntDigits) else // большое целое число. Определяем экспериментально /// вариант SignificantPower := 17 - IntDigits N:=abs(N) SignificantPower := -16 while SignificantPower <= 16 if !(N == N+10^(SignificantPower)/2) exit end ++SignificantPower end end ret SignificantPower *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func FloatPrecision(N) // определение "маленького числа" по заданной во FloatDigits() длине разрядной сетки local Precision:=0 //ret FloatAddIsSensitive(N) // это эквивалентный вариант (раньше при округлении до дробных разрядов. проверить на целых разрядахRoundTo(...,+1) Precision:=10^(SignifPower(N))/2 ///new!!! в 2 раза меньше! ret Precision *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func FloatAddIsSensitive(N) // экспериментальное определение, прибавление какого маленького числа начинает быть чувствительным для данного числа local i:=-16,Sign Sign := if(N>=0,1,-1) N:=abs(N) while i <= 16 if !(N == N+10^(i)) ret 10^(i) * Sign * 1 end ++i end ret 0 *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func IsInt(N) ret abs(N-int(N)) <= FloatPrecision(N) *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func IsEven(N) ret abs(N % 2) <= FloatPrecision(N) *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func RealEQ(N1,N2) ret FloatEQ(N1,N2) func FloatEQ(N1,N2) ret abs(N1-N2) <= FloatPrecision(N1) *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func RoundTo(N,Dec,Round5toUp) local R,Sign,RoundDigit,PrevDigit,DigitIsLast Default(Round5toUp,Round5DefaultMode) Sign := if(N >= 0,1,-1) //N += Sign*FloatPrecision(N) RoundDigit:= N RoundDigit:= RoundDigit*10^(Dec+1) // округляемая цифра сдвинута на позицию единиц /* if kbState(ScrollState)#0 xview({ RoundDigit,; int(RoundDigit),; abs(RoundDigit-int(RoundDigit)-1),; FloatPrecision(RoundDigit)}) end */ if Dec >=1 .and. abs(RoundDigit-int(RoundDigit)-1) < FloatPrecision(RoundDigit) // нельзя FloatEQ(RoundDigit-int(RoundDigit),1) т.к. точность сравнения будет взята, ориентируясь на результат операции (единицу), а не на исходные числа // После запятой одни девятки // N==0.5035==0.5034(999) выполняем RoundTo(0.5035,3) // RoundDigit - почти целое число // N*10000=5034.999999999999, int(N*10000)=5034, а не 5035 как должны было бы быть // N-int(N) == 0.9999999999990905 почти единица // ?"!!!!!",str(N,60,30),CRLF // __wait("") DigitIsLast := .T. while abs(RoundDigit-int(RoundDigit)-1) < FloatPrecision(RoundDigit) RoundDigit+=FloatPrecision(RoundDigit) nop end else // дробная часть оч маленькая DigitIsLast := (abs(RoundDigit-int(RoundDigit)) < FloatPrecision(RoundDigit)) end RoundDigit:= int(RoundDigit) RoundDigit:= val(substr(nTrim(RoundDigit),-1)) if RoundDigit == 5 .and. DigitIsLast PrevDigit := int(N*10^Dec) PrevDigit := val(substr(nTrim(PrevDigit),-1)) if !isEven(PrevDigit) .or. Round5toUp N += 10^(-1*(Dec+1)) * Sign else N -= 10^(-1*(Dec+1)) * Sign end end R := round(N,Dec) ret R *▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄▄ *█ █ *▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀▀ func ScienceFormat(N) local i,IntDigits,Power:=0,S,Sign Sign := if(N<0,.T.,.F.) N := abs(N) if N < 1 while N < 1 Power-- N*=10 end else IntDigits := len(nTrim(int(N))) Power := IntDigits-1 N := N/(10^Power) end S := if(Sign,"-","")+str(N,3,1)+"E"+nTrim(Power) ret S [/pre2] SUV.CH [pre2] #define nTrim(N) ltrim(str(N)) #define Default(X,V) if(X == NIL,X := V,) #command SetGet Func <F> => func <F>(X) ;; local Old ;; Old := <F> ;; if X # NIL ;; <F> := X ;; end ;; ret Old #command SetGet Func <F> for <V> => func <F>(X);; local Old ;; Old := <V> ;; if X # NIL ;; <V> := X ;; end ;; ret Old #xcommand ret => return #xcommand ret <X> => return <X> #command repeat => while .T. #command until <ExitExpr> => if <ExitExpr> ; exit ; end ; enddo #command nop => if .T. ; end [/pre2]
NickSam: Dima, спасибо большое! Попробовал P_Round и Ft_round, у меня проблема не ушла. Почитал комментарии к RoundTo и кажется понял в чем дело. Для вещественных чисел Harbour отводит 17 знаков для хранения (не только под дробную часть, но и под целую). Вот пример где надо округлять до 11 знаков после запятой: set decimals to 19 m_cena = 3537.65 / ( 1 + (10/100) ) alert(m_cena) alert(P_Round(m_cena,11)) Результат: 3216.0454545454550000000 3216.04545454546 Результат этого же вычисления на Python: 3216.0454545454545041139 3216.04545454545 Можно ли как то побороть это на Harbour? Может есть какая то функция на C где результат хранится с большей точностью?
Pasha: Харбор использует 8-ми байтовый тип double, в котором под мантиссу отводится 52 бита 2^52 получаем 4503599627370496, т.е. получаем точность не более 16 десятичных знаков. Говорить о какой-то точности и правильности округления сверх этого бессмысленно, это ограничение типа double Откуда питон берет значность сверх этого ? Есть еще 16-байтовый формат чисел с плавающей точкой, т.е 128 бит, может он и используется. Но харбор использует 64-битный формат
Pasha: Касательно расширения стандартного типа double 64-бит (8 байт), который поддерживается в Харборе, с мантиссой 52 бита, экспонента 11 бит, и знаковым разрядом (52+11+1=64). Точность 16 десятичных знаков, как мы видели. Это суммарные десятичные знаки, перед запятой плюс после запятой. 128-битый формат, который я упоминал (16 байт) - у него мантисса 112 бит, степень 15 бит, и знаковый разряд В C компиляторах он обозначается как long double. Но везде по разному В некоторых компиляторах так: bcc32: long double это 10 байт, или 80 бит. Мантисса 64 бита, степень 15 бит, и знак. Это дает точность 19-20 десятичных знаков: 2^64 = 18446744073709551615 mingw32: здесь long double это 12 байт, или 96 бит. Сколько там мантисса не нашел, но если порядок степени везде 15, то мантисса должна быть 80 бит. Это дает точность 24 десятичных знака В mingw64 long double это как раз 128 бит, или 16 байт. Мантисса 112 бит, степень опять 15 бит, и знак. 2^112 дает не знаю сколько точно, такой архитектуры еще нет, чтобы посчитать. Надо делать отдельную программу для вычисления 2^112. Но там точность 33-34 десятичных знака. Но все эти расширения вычислений с плавающей точкой сверх 64 бит аппаратно не поддерживаются, это программная реализация. Расчеты намного медленнее. Хотя с современной производительностью процессоров это будет незаметно В связи с этим мне вспоминаются времена, когда сопроцессор это было отдельное устройство, и если его не было, была только программная эмуляция, и это тогда было заметно, и очень 2-3 года назад мне пришлось конвертировать базу из очень старой программы на Турбо Паскале, там был такой интересный тип данных 6 байт real, это 48 бит. Разбирал его побитно. По сравнению с типом float 32 бита это был прогресс. Одно непонятно: для чего может понадобиться такая точность. Даже госдолг США (31 терадоллар) пока вписывается - надо 14 знаков. Правда если считать центы, то 16, это уже на пределе, скоро точности double будет не хватать. Сделал программу, которая эмулирует long double через строки харбора, и базовые арифметические операции для этих чисел. Результат кстати получается не такой, как дает питон. Похоже питон тоже дает точность порядка 16 разрядов, а сверх того - это цифровой шум Ссылка: https://cloud.mail.ru/public/V6V4/VPHY8RNHA
SergKis: Pasha У нас это не доступно PS Когда то все считали в целых числах, выводили как строки
Pasha: Нет доступа к mail.ru ? Текст не очень большой, брошу здесь. [pre]proc main Local cld1 := ld_set( 3537.65 ) Local cld2 := ld_set( 1.1 ) Local cld := ld_div( cld1, cld2 ) ? len(cld1) ? ld_get(cld1) ? ld_get(cld) ? ld_str(cld2, 8) ? ld_str(cld,, 18) ? ld_str(cld, 24, 18) ? ld_str(cld, 40, 34) retu #pragma BEGINDUMP #include "hbapi.h" #define ld_len sizeof( long double ) typedef union { long double value; char string[ ld_len ]; } xld16; HB_FUNC( LD_SET ) { xld16 ld; ld.value = hb_parnd( 1 ); hb_retclen( ( char * ) &ld, ld_len ); } static xld16 ld_get( int iPar ) { HB_SIZE ulLen = hb_parclen( iPar ); const char *szLD = hb_parc( iPar ); xld16 ld; if( ulLen == ld_len ) memcpy( &ld.string, szLD, ld_len ); else ld.value = 0.0; return ld; } HB_FUNC( LD_GET ) { xld16 ld = ld_get( 1 ); hb_retnd( ld.value ); } HB_FUNC( LD_ADD ) { xld16 ld1 = ld_get( 1 ); xld16 ld2 = ld_get( 2 ); xld16 ld; ld.value = ld1.value + ld2.value; hb_retclen( ( char * ) &ld, ld_len ); } HB_FUNC( LD_SUB ) { xld16 ld1 = ld_get( 1 ); xld16 ld2 = ld_get( 2 ); xld16 ld; ld.value = ld1.value - ld2.value; hb_retclen( ( char * ) &ld, ld_len ); } HB_FUNC( LD_MUL ) { xld16 ld1 = ld_get( 1 ); xld16 ld2 = ld_get( 2 ); xld16 ld; ld.value = ld1.value * ld2.value; hb_retclen( ( char * ) &ld, ld_len ); } HB_FUNC( LD_DIV ) { xld16 ld1 = ld_get( 1 ); xld16 ld2 = ld_get( 2 ); xld16 ld; ld.value = ld1.value / ld2.value; hb_retclen( ( char * ) &ld, ld_len ); } HB_FUNC( LD_STR ) { int iLen = hb_parni( 2 ), iDec = hb_parni( 3 ), iLen1; xld16 ld = ld_get( 1 ); char szStr[50]; char sFormat[10]; int iFmtLen = 1; sFormat[0] = '%'; if( iDec > 0 ) { char sDec[10]; sFormat[1] = '.'; itoa( iDec, sDec, 10 ); memcpy( sFormat + 2, sDec, strlen( sDec ) ); iFmtLen += 1 + strlen( sDec ); } sFormat[ iFmtLen ] = 'L'; sFormat[ iFmtLen + 1 ] = 'f'; sFormat[ iFmtLen + 2 ] = 0; #if defined(__MINGW32__) || defined(__MINGW64__) __mingw_snprintf(szStr, 50, sFormat, ld.value ); #else snprintf(szStr, 50, sFormat, ld.value ); #endif iLen1 = strlen( szStr ); if( iLen > 0 && iLen > iLen1 ) { if( iLen1 < iLen && iLen < 39 ) { int i; for( i = iLen1 - 1; i >= 0; i -- ) szStr[ i + ( iLen - iLen1 )] = szStr[ i ]; memset( szStr, ' ', iLen - iLen1 ); iLen1 = iLen; } } hb_retclen( szStr, iLen1 ); } #pragma ENDDUMP [/pre]
SergKis: Pasha пишет Нет доступа к mail.ru ? И не только, ко многому заблокировано. С впн в облако попадаешь (оооочень долго), но взять не дает Всем, кто имел почту на mail.ru, пришлось делать новые адреса на доступных ресурсах брошу здесь Оч. похоже как операции с целыми, точку ставили при выводе в нужную позицию
Pasha: Можно оформить как класс с перегрузкой операторов + - * /, тогда использование этого типа будет выглядеть почти как работа с обычными числами
полная версия страницы